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PREFACE 


The first experimental studies of laser-matter interaction were conducted in 
the early 1960s when the first high-power lasers were created. In these works, 
mainly thermal effects were studied including heating, melting, evaporation 
of condensed matter, and thermionic electron emission. When higher laser 
intensities were reached, it became possible to study the optical breakdown 
of gas and solid dielectrics. In the 1960s, the first theoretical models of laser- 
matter interaction also were proposed. The simplest models were related to 
moderate laser intensity ranges (I ~ 10 5 -10 9 W/cm 2 ) and took into account 
the energy transfer in condensed matter, the kinetics of phase transitions, 
and the dynamics of evaporated material expansion. 1-4 A survey of the re¬ 
sults obtained during the initial period of laser-matter interaction studies can 
be found in Ready 5 and Anisimov et al. 6 Despite many simplifications, the¬ 
oretical models 1-4 provide a basic outline of the actual interaction of laser 
radiation with solids at moderate laser intensities. However, it appeared to be 
difficult to reach a quantitative agreement between theory and experiments. 
In particular, the detailed analysis of phase composition and energy balance in 
laser ablation of metals 1,5,7 shows that a considerable part of the ablated ma¬ 
terial is in the form of liquid droplets, and the specific ablation energy—equal 
to the ratio of the incident laser energy to the mass of the material removed 
from the target—is several times less than the specific heat of vaporization. 
These facts cannot be explained by simple theoretical models. Another mani¬ 
festation of the very complicated character of actual laser-matter interaction 
processes and the inadequacies of simple theoretical models follows from stud¬ 
ies of target reflectivity changes during laser irradiation. 5,8 The experiments 
show that total reflectivity and particularly its specular part, is reduced con¬ 
siderably during a laser pulse. This reduction cannot be explained easily by 
a decrease in surface layer conductivity due to heating. Finally, such im¬ 
portant characteristics as the phase composition of the ablated material, the 
specific ablation energy, and the recoil momentum transferred to the target 
are extremely sensitive to the specifics of the laser pulse structure and focusing 
conditions. 9 

These facts can be understood from a single viewpoint if we suppose that 
laser-induced evaporation of condensed materials can be unstable under cer¬ 
tain conditions. The first example of laser-produced evaporation instability, 
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studied by Anisimov, Tribel’skii, and Epel’baum 10 was the corrugation insta¬ 
bility of a plane stationary vaporization front in highly absorbing solids. This 
instability is driven by the temperature gradient in the vicinity of the phase- 
transition front. A similar instability mode appears in laser-induced evap¬ 
oration of transparent dielectrics. Furthermore, the thermal and avalanche 
breakdown of transparent dielectrics also can be considered to be a partic¬ 
ular mode of instability of these materials in a strong electromagnetic field. 
Many other examples of instability have been studied during recent years. 11-13 
These studies clearly demonstrate that most of the laser-matter interaction 
processes are unstable. 

In the present book, a review is given of thermal and hydrodynamic insta¬ 
bilities appearing in laser-matter interactions at moderate intensities. These 
instabilities also can arise in the processes of heating of condensed matter 
with electron and ion beams, shock waves, etc. Nonresonant interactions are 
the subject of this book; to include resonant electromagnetic modes related to 
the generation of surface electromagnetic waves would require a second book. 
Note also that these phenomena are not characteristic of the interaction of 
electron or ion beams with solids. A comprehensive review of the phenomena 
associated with surface electromagnetic waves can be found in references 11- 
13. We will not consider here problems of laser-driven inertial confinement 
fusion; instead, we will restrict our attention to the range of low and moderate 
laser intensities that are important for technological applications of lasers. 

This book is constructed as follows: In the first four chapters, we give a 
survey of the basic processes of nonresonant laser-matter interactions. Laser- 
induced breakdown of transparent dielectrics is considered in Chapter 2, as an 
example of instability of “thermal explosion” type (see Frank-Kamenetskii 14 ). 
Chapter 4 contains a brief analysis of the effects produced by ultrashort (pi¬ 
cosecond and femtosecond) laser pulses. This problem recently has attracted 
extensive interest due to the important applications of ultrafast laser tech¬ 
nology. The problems of stability of thermal processes induced by ultrashort 
laser pulses have not yet been studied. 

In chapters 5 to 8 we consider thermal and hydrodynamic instabilities. 
This provides a theoretical background to the interpretation of experimental 
results and an understanding of the effect of instabilities on the processes of 
laser technology. 

We would like to note that the instabilities discussed in this book are 
important for the interpretation of laser interaction experiments and techno¬ 
logical applications of lasers. The instabilities are also of interest from the 
standpoint of fundamental physics. There are many examples of physical sys¬ 
tems whose stable steady states have a symmetry that does not correspond to 
the symmetry of external conditions. In these cases, it is customary to discuss 
spontaneous symmetry breaking. Many examples of this general phenomenon 
are known in the fields of atomic and molecular physics, solid state physics, 
hydrodynamics, and quantum field theory. 15-17 Laser-matter interactions give 
additional examples of symmetry breaking due to instability. The analysis of 
the symmetry breaking and structure formation is a part of the general theory 
of self-organization, which has been given notice during the past 20 years. The 
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concepts and methods of this theory can be employed, as we will see later, in 
laser-matter interaction studies. 

This book is devoted mainly to the theoretical analysis of laser-matter 
interaction and instability phenomena. This analysis is based primarily on 
the results of studies carried out by Russian physicists during the past 10-15 
years. These results were published primarily in Russian journals and were 
delayed in reaching western readers. 

The authors are grateful to their friends and colleagues: D. Bauerle, 
A. M. Bonch-Bruevich, M. N. Libenson, B. S. Luk’yanchuk, J. Meyer-ter- 
Vehn, M. I. Tribel’skii, and S. Witkowski for helpful discussions and criticism 
during the preparation of this book. One of the authors (S. A.) acknowledges 
with gratitude the kind hospitality of the Max-Planck-Institute for Quantum 
Optics (Garching, Germany) and the support of the Alexander von Humboldt 
Foundation. 
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Chapter 1 

INTERACTIONS 
OF LASER RADIATION 
WITH HIGHLY ABSORBING 
MATERIALS 


Generally, the entire field of laser-matter interactions can be divided into two 
parts: (1) resonant and (2) nonresonant interactions. The first part includes 
selective excitation processes in systems with discrete spectra, multiphoton 
ionization of atoms and molecules, multiphoton photoelectron emission, and 
nonlinear optical phenomena in solids and liquids. Such factors as coher¬ 
ence, polarization, and spectral characteristics of laser radiation affect these 
excitation processes. Laser-induced thermo chemical and thermophysical phe¬ 
nomena can be related to nonresonant interactions, as well as laser breakdown 
of gases and solids and laser plasma generation. For this group of phenomena, 
the coherence and spectral characteristics of the laser radiation are not very 
important. Principal effects are connected, in this case, with the heating of 
a material, and the important factors are the laser intensity, pulse duration, 
and focusing conditions. In this book, we shall discuss mainly the thermal 
effects of laser radiation. The emphasis will be on laser heating, melting, and 
vaporization of solid materials and the related instability of these processes. 
We consider also the dynamics of vapor plume expansion, with an emphasis 
on hydrodynamic instabilities of vapor expansion. 


1.1 Laser heating of highly absorbing solids 
without phase change 

When laser radiation intensity is rather low, no phase transition occurs, and 
the only effect of laser absorption is heating of the material. We first consider 
the laser heating of metals, the most important materials with high optical 
absorption. In metals, the laser radiation is absorbed by “free” electrons, 
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and energy transfer in metals is also due to electron heat conduction. The 
temperature field is described by the standard heat conduction equation: 

pcdT/dt = div (/cVT) + Q{ r, t) (1.1) 

where T is the temperature, k is the thermal conductivity, c is the specific 
heat, p is the density, and Q(r,£) is the heat production (per unit volume per 
unit time) due to the laser radiation absorption. In the general case, c and 
k are temperature-dependent, and Equation (1.1) is nonlinear. For metals, 
however, at temperatures above the Debye temperature, c and k can be taken 
as constants. 1,2 

It should be noted that the macroscopic approach to the problem of laser 
heating is justified when the characteristic time and space scales of the tem¬ 
perature field are much larger than the mean free time and mean free path of 
electrons, respectively. These conditions are fulfilled in the case of nanosecond 
(and longer) laser pulses. However, in the case of picosecond and subpicosec¬ 
ond pulses, the above approach should be modified. 

The heat source Q(r,t) in Equation (1.1) depends on the laser pulse pa¬ 
rameters and on the optical properties of the material irradiated. In many 
practical cases, the Drude-Zener theory 2,3 provides an adequate description of 
the optical properties of materials whose optical absorption is mainly due to 
free electrons. Solid and liquid metals are typical examples of such materials. 
The Drude-Zener theory leads to the following expression for Q(r,£): 

Q = Apl(x , y, t) exp(-/iz) (1.2) 

where A is the surface absorptivity, /i is the absorption coefficient of the 
material, and /(#, y, t) is laser radiation intensity at the material surface (z = 
0). The absorption coefficient, /i, is related to the skin depth, <5, by /i = 2/6. 
According to Sokolov 2 and Libenson et al. 4 /i is independent of temperature, 
while A — 1 — R (R is the surface reflectivity) is a linear function of the surface 
temperature: 

A(T) = Ao + Ai(T — To) (1.3) 

where Ao is the surface absorptivity at room temperature, To- This tem¬ 
perature dependence of A results from the fact that A is proportional to the 
electron-phonon collision frequency which, in turn, is proportional to the crys¬ 
tal lattice temperature. 

Boundary conditions for Equation (1.1) can be written in the following 
form. At a large distance from the material surface, as z —> oo 

T —> To = constant (1.4) 

(in the majority of cases, To can be set equal to zero). On the material surface, 
z = 0, the continuity condition for the energy flux should be imposed. It reads 
as: 

—ndT/dz = q c + q r 

where q c and q r are energy fluxes due to convection and radiation, respectively. 
In the case under consideration (no phase transition) the temperature of the 
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surface is not very high, and the energy losses are negligible, compared to the 
laser energy flux. An estimate of the region of laser-pulse parameters in which 
the energy losses are negligible is given by Ready . 1 We can thus set: 

dTjdz — 0 at z = 0 ( 1 . 5 ) 

The initial condition for Equation ( 1 . 1 ) reads: 

T(r, 0 ) = T 0 ( 1 . 6 ) 

In many cases of practical interest, the transverse dimensions of the laser 
focusing spot are large compared to the thickness of the heated layer, and 
the heat conduction problem (Equations ( 1 . 1 ) and (1.4) through ( 1 . 6 )) can be 
considered one-dimensional. 

Note that the heat conduction problem (1.1) with Q( r, t) defined by Equa¬ 
tions (1.2) and (1.3) is linear and can be solved using standard methods . 5 
In the simplest case of time- and coordinate-independent laser intensity, 
/(x, y, t) = Jo = constant, this problem has been solved by Libenson et al . 4 
To obtain the solution, we apply the Laplace transform, defined as: 


OG 


u(z,p) = J T{Zi t) exp(— pt) dt 


o 


to Equation (1.1) and to the boundary conditions (1.4) and (1.5). The ordi¬ 
nary differential equation that results from the Laplace transformation can be 
solved readily to obtain: 


u(z,p) = (plo/cp) [(Aiuop 4- A 0 )/p(p - XP 2 )] 

x [exp(—/iz) - (p/P) exp(—/3z)] 


(1.7) 


Here, uo = u(0,p), P = (p/x) 1 ^ 2 ^ an< 3 X — K / C P ls the heat diffusivity. The 
temperature, T(z,£), is calculated by applying the inverse Laplace transform 
to Equation (1.7). To avoid cumbersome formulas, we perform the calculation 
only for the surface temperature, T( 0 ,£). Solving Equation (1.7) for uo, we 
obtain 


u 0 = [J 0 A 0 /i/pcp(ai - a 2 )] 


_ 1 _ 1 

y/p-Oil y/p-OL* 


(1.8) 


Here ai >2 = (Py/x/ 2)(—1 ± y/1 -f 5), and S = 4/oAi//i/t. The parameter 
S is proportional to the laser intensity, and in the case under consideration, 
it is relatively small. For typical metals (A\ ~ 10~ 4 K" 1 , fi ~ 10 5 cm -1 , 
k = 3 W/cm K), S ~ 10~ 2 at laser intensity 10 7 W/cm 2 . Assuming S < 1 
and performing the inverse Laplace transform, we obtain from (1.8): 

T(0,t) = To + (IoAq/K fi) [exp(/x 2 x0 erfc(fi^/xt) — l] 9 ) 

+ (Ao/Ax) [exp( 7 t)erfc(- v / 7 ?) - 1] 


OO 


erfc(x) = (2/y/7r) J exp(—y 2 )dy, and 7 = IqA 2 /kcp 


where 
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In the majority of calculations regarding the laser heating of solid surfaces, 
the temperature dependence of the surface reflectivity usually is neglected. 
This results in the following equation: 

T( 0, f) = T 0 + (IqA 0 /Kfi) [2 Hy/xt/y/ir + exp(/x 2 xt) erfc(^v^) - l] (1.10) 

which follows from (1.9) if we set A\ = 0 . It is easy to see that Equations 
(1.9) and (1.10) give almost the same values of the surface temperature if 
jt -C 1. However, when t 7"* 1 , the changes in surface reflectivity become 
important. This leads to exponential time growth of the surface temperature 
at large t: 

T(0,t) ~ To 4- 2(j4oMi)exp(7t) (1.11) 

Note that when 5 1 , only the regime with exponential temperature growth 

exists, described by the following equation: 

T( 0 ,t) = (Aq/A\) [exp(I 0 Aifxt/c) - 1 ] 

In many cases of practical interest, the skin depth is much smaller than the 
spatial scale of the temperature field, and we can consider light absorption as 
a surface process. Thus, we can set Q = 0 in ( 1 . 1 ) and change the boundary 
condition (1.5) at z = 0 to: 


— ndT/dz = A(T)Iq ( 1 . 12 ) 

The solution of the last problem can be obtained formally from (1.9) by setting 
fi —> oo. The result is given by the equation: 

T(Q,t) = (A 0 /Ai) [exp( 7 t) erfc(—v/ 7 t) - l] 

which coincides with ( 1 . 11 ) in the limiting case 7 1 1 . 

The spatial temperature profile can be calculated by taking the inverse 
Laplace transform of (1.7). Since the result for a general case is cumbersome, 
we will consider only the simplest case of constant reflectivity, A\ = 0 . The 
temperature field is given by : 5 

T(z, t) = (2loAo/K)y/xt ieiic(z/2y/xt) - (IqAq/hk) exp(-/xz) 

+ (IoAo/2fm) exp(/i 2 xt - V>z) erfc(/i x /xt - z/2y/xt) ^ 13 ) 
+ (/ 0 j4o/2/iac) exp(/i 2 xt + l*z) erfc (tiy/xt + z/2y/xt) 

where 

00 

ierfc(x) = J eiic(y)dy. For z = 0 , (1.13) reduces to ( 1 . 10 ). 

X 

As fi tends to infinity (surface absorption), we have from (1.13) the following 
simple result: 

T(z,t) = (2J 0 i4o/K)v^ierfc(z/2Vx?) 

We have considered the one-dimensional temperature field produced by 
laser heating of a highly absorbing solid. Note that this field is nonstationary: 
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the surface temperature increases over a long time as exp( 7 ^) (at S 1), 
or as exp( x /x 7 fit) (at S 1), and the thickness of the heated surface layer 
increases as y/xt. 

A one-dimensional approximation is valid when the heated layer thickness 
h (x y/xt is much smaller than the laser spot size R. It is clear that this 
condition is broken when t > R 2 /x • Thus, when the laser pulse length is larger 
than R 2 /Xi the temperature distribution becomes two- or three-dimensional. 
This situation occurs, for example, when a material is heated by a continuous 
laser. In this case, the spatial scale of temperature distribution is much larger 
than the skin depth, and the light absorption may be considered as a surface 
process. We can, therefore, set Q = 0 in Equation (1.1) and use the boundary 
condition similar to (1.12) with laser flux intensity Jo, localized within the 
laser focusing spot. It should be noted that energy loss due to convection and 
radiation can be no longer neglected, since for long laser pulses, a large area 
around the laser spot becomes hot and contributes to heat loss. In this case, 
the heat loss from the surface z = 0 may reach the same order of magnitude as 
the laser power absorbed. The linear heat conduction equation (1.1) with Q = 
0 and linear boundary condition at z = 0, can be solved easily using standard 
methods. 5 Examples of multidimensional temperature field calculations can 
be found in Ready 1 and Carslaw and Jaeger. 5 We will not consider the details 
of these calculations; note that for a time-independent laser intensity a steady- 
state temperature profile can be reached in the multidimensional case. 

The above analysis of heat transfer in a laser-heated solid contains a num¬ 
ber of simplifications. We did not consider the real temporal shape of a laser 
pulse, which can be important for quantitative calculation of the tempera¬ 
ture. In the framework of linear heat conduction problems, the consideration 
of arbitrary laser pulse shape does not lead to any complications, but gives 
rise to more laborious calculations, which sometimes require the application 
of numerical methods. More serious difficulties may be related to the nonlin¬ 
earity of the heat transfer problem that results from temperature dependence 
of thermophysical and optical characteristics of the material. Nonlinearity 
also can appear in the boundary condition if radiative energy loss is taken 
into account. Nonlinear heat conduction problems may be solved, generally, 
by approximate or numerical methods. Examples of such calculations can be 
found in Prokhorov et al. 6 


1.2 Laser-induced melting 

We now turn to higher laser energy absorption and consider the melting of 
solids absorbing laser radiation. In practice, laser-induced melting usually is 
accompanied by vaporization. Since the vaporization rate depends strongly 
on the temperature, melting without substantial vaporization can be observed 
only in a very narrow range of laser parameters. As noted by A. I. Shal’nikov, 
this range depends on the saturated vapor pressure of a particular material 
at its melting point. 7 It is important to note that for different materials, this 
parameter varies by many orders of magnitude. For example, the saturated 
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vapor pressure at the melting point equals 9.8 x 10" 3 bar for Cr, 5.7 x 10" 2 bar 
for Br, and 2.4 x 10“ 1 bar for Gd; at the same time, for Ga this pressure is 
9.2 x 10" 41 bar, for Sn 5.7 x 10“ 26 bar and for S 2.6 x 10“ 25 bar. 8 It is clear that 
for the materials with rather high saturated vapor pressures at the melting 
point, pure melting without vaporization rarely is observed. 

The simplest estimate of the mass of fused material can be obtained if 
one supposes that the melting front coincides with the melting temperature 
isotherm. In this estimate, the latent heat of fusion is neglected, which is not 
small, in most cases, in comparison with the total enthalpy of a material at 
the melting point. A more accurate calculation of melting front position can 
be made in terms of the so-called “problem of Stefan.” 5 With this approach, 
the temperature fields in the solid and liquid phases are described by the 
heat conduction equations, and the boundary condition on the moving phase 
boundary is assumed to have the following form: 

- KidTi/dn + n s dT s /dn = V n pLm (1.14) 

where d/dn is the normal derivative, L m is the latent heat of fusion, and V n is 
the melting front velocity. The melting front is defined, in the framework of the 
“problem of Stefan,” as an isothermal surface where the temperature is equal 
to the melting temperature, T m . Note that the law of motion of the melting 
front is not known a priori and should be found as an eigenvalue in the course 
of solution of the heat transfer problem. An important example of such a 
calculation is given by Carslaw and Jaeger. 5 This is the well known Neumann’s 
solution for the semi-infinite region, initially at a constant temperature To > 
T m . Physically, this solution corresponds to instantaneous heating of the 
solid material and may be considered, in principle, as a model for pulsed-laser 
produced melting. It should be noted, however, that Neumann’s solution 
leads to a physically meaningless conclusion that the melting front velocity 
becomes infinite, V n = oo, at finite superheating, (To — T m )c/L m > 1. This 
qualitatively incorrect result is the consequence of the basic assumption of 
the “problem of Stefan” that the melting front dynamics are governed by 
the energy flux only and are independent of the phase transition kinetics. It 
is clear that this simple assumption is applicable only if the melting front 
velocity is very small. 

Melting kinetics at low superheating have been considered by Motorin 
and Musher. 9 Assuming the melting front velocity to be proportional to the 
difference of chemical potentials of solid and liquid phases, fi s — m, Motorin 
and Musher 9 obtained the following formula for V n : 

Vn = [c(T - T m )/L m + (cTmc/pL^il/R! + l/Ri)] (1.15) 

where a is the surface energy density (surface tension) at the solid-liquid 
interface, R\ and R 2 are the principal radii of curvature of the phase boundary, 
and V m is an empirical constant that can be determined experimentally. In 
order of magnitude, V m is close to the speed of sound. For copper, e.g., V m 
equals 7 x 10 4 cm/s. 9 When phase transition kinetics are taken into account, 
the temperature at the melting front can no longer be assumed to be equal 
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to the normal melting point. The solid phase at the front has a temperature 
higher than that of the melting point, and the front velocity, as seen from 
(1.15), is proportional to the superheating. To calculate the melting front 
velocity and temperature, one should solve the heat conduction equations for 
the solid and liquid phases with boundary conditions (1.14) and (1.15) imposed 
on the front. These calculations (for a more general case with melting and 
evaporation kinetics both taken into account) were performed by Anisimov 
and Barsukov. 10 We shall discuss these calculations in Section 1.4. 

It is important to note that the melting can be considered to be a surface 
process only when the superheating is sufficiently low. In the general case, the 
bulk formation of liquid phase nuclei (homogeneous nucleation) must be taken 
into account. There is a fundamental distinction between bulk and surface 
melting. 11 Usually, a liquid completely wets the surface of a solid phase of the 
same substance. Therefore, the formation of a liquid nucleus on the surface 
does not require work to be done to form a new surface. In other words, there 
is no energy barrier for the surface melting. On the contrary, the formation of 
a liquid nucleus within a crystal involves elastic deformations; this makes the 
bulk nucleation thermodynamically disadvantageous. This is true, of course, 
for any first-order phase transition in solids. The existence of the finite energy 
barrier for the bulk nucleation is used sometimes to explain the fact that first- 
order phase transitions do not begin on the phase equilibrium line. 

It should be noted that the work required to form a nucleus of a new phase 
depends on the nucleus shape. For example, in the case of melting, a spherical 
nucleus formation can be expected to occur at the superheating of the order: 12 

_ 2 cT m KG(3 2 

(4G 4- ZK)p s L 2 m 

where K is the compressibility of the liquid phase; G is the shear modulus of 
the solid phase; and f3 = ( p s — pi)/ p 5 , p 5 , pi is the density of solid and liquid 
phase, respectively. It was shown by Motorin and Musher, 12,13 however, that 
for disk-shaped nuclei, the effect of elastic deformations on the probability of 
nuclei formation is negligible. Thus, bulk melting can be expected to occur 
when the superheating is rather small. The value of superheating here is 
related to the finite growth time of the nucleus. 


1.3 Laser-induced evaporation 

As mentioned in the preceding section, melting without substantial vaporiza¬ 
tion can be produced only in a very narrow range of laser parameters. Usually, 
both phase transitions occur simultaneously under conditions typical for laser 
experiments. Note that the latent heat of vaporization is much larger than 
that of fusion (typically, by a factor of 20-50 times). Thus, evaporation plays 
the most important part in the energy balance. 

In this book, we examine the interactions of laser radiations with matter 
at moderate laser intensities, I ~ 10 5 -10 8 W/cm 2 . At these intensities, the 
temperature is usually well below the thermodynamic critical point of the ma¬ 
terial. Under this condition, a sharp boundary exists between the condensed 



16 


INSTABILITIES IN LASER-MATTER INTERACTION 


and gaseous phases, its thickness being on the order of interatomic distance, 
and the density of the vapor is several orders of magnitude less than in the 
condensed phase. In the majority of cases of practical interest, light absorp¬ 
tion in the vapor is relatively small. To simplify our discussion we will begin 
with the assumption that the evaporation occurs in vacuum and the vapor 
does not absorb the laser radiation. Note that for some materials, the last 
assumption is not quite correct, especially near the upper boundary of the 
previous laser intensity range. 

When the vapor plume does not absorb the laser light, the processes in the 
condensed and gaseous phases can be considered separately. As in the previous 
chapter, we shall describe the temperature field in the condensed phase by the 
heat conduction equation (1.1) with boundary conditions taking into account 
phase transition kinetics and energy balance. Note that the vaporization of 
condensed matter cannot be described, even in a very rough approximation, 
in terms of the “problem of Stefan” with a fixed phase transition temperature. 
Actually, a solid (liquid) evaporates at any temperature not equal to 0 K. The 
evaporation rate strongly depends on the temperature. This dependence, in 
the case under consideration, can be written in the form: 14,15 

V n = V 0 exp (-U/T) (1.16) 

where V n is the normal component of the evaporation front velocity, U — 
MLyjks , L v is the latent heat of vaporization (per unit mass), ks is Boltz¬ 
mann’s constant, M is the atomic mass, and Vo is a constant whose value is 
of the order of the speed of sound in the condensed phase. We will show later 
how the evaporation rate is related to the saturated vapor pressure. Detailed 
data on saturated vapor pressures and evaporation rate of different materials 
can be found in references 8, 16, and 17. 

Now we consider the simple problem of plane evaporation wave propaga¬ 
tion. 14,18 Let the boundary between the condensed and gaseous phases be the 
plane and its law of motion due to condensed phase vaporization be z = Z(t). 
Introducing new independent variables: 

C — ii[z — Z(t)] and r — fi 2 xt 
we write the heat conduction equation (1.1) in the form 

dT/dr = d 2 T/d{ 2 + /3(r)dT/dC + Q( C, r) (1.17) 

where f3(r) = V //ix; V = dZjdt\ Q((,t) = Kexp(—(); and K = I(t)A/fiK. 
The boundary conditions on the moving evaporation front £ = 0 are: 

dT/d C = (A H/c)(3(t) and (3(r) = (V 0 /fxn) exp[-[//T(0, r)] (1.18) 

Here, AH is the difference in enthalpies of the gaseous and condensed phases, 
with respect to unit mass. It is approximately equal to AH ~ L v +ksT/2M ~ 
Ly . 

If the laser intensity I is constant, a steady-state solution of the boundary 
value problem (1.17) and (1.18) exists. It describes the stationary temperature 
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Figure 1.1: Temperature and velocity of the stationary evaporation fronts as 
functions of laser intensity. 


profile T S (Q propagating with constant velocity V s through the solid (liquid) 
material. Assuming in (1.17) and (1.18) that /3(r) = (3 = constant; T(£,t) = 
T S (Q; and dT s /dr = 0, we arrive at an ordinary differential equation, which 
can be solved as follows: 

T s ( C) = exp(-C) 4- C 2 exp(-^C) (1.19) 

where C\ = K/(f3 — 1) and C 2 = —AH/c — K//3(/3 — 1). It follows from (1.19) 
that the derivative, dT s /d£ at £ = 0, is positive and equal to L v V s /ck. Since 
at £ —* oo, T 5 (C) —* 0, and dT s /d£ < 0, we conclude that the steady-state 
temperature at the phase boundary is lower than in the bulk of material. This 
fact plays a decisive role in the stability analysis, which we shall describe in 
Chapter 5. 

If we set C = 0 in (1.19) and use the equations for Ci, C 2 , and K, we 
obtain the energy conservation equation, which can be written as: 

I abs — AI — V s p[L v 4 cT s (0)] (1.20) 

where I abs is the absorbed laser intensity. Taking the equation of vaporization 
kinetics (1.15) in the form: 


V s =V 0 exp[-U/T s (0)} (1.21) 

and solving the set of equations (1.20) and (1.21), we determine the stationary 
velocity, V s , and the stationary temperature, T s ( 0), on the surface of the 
condensed phase. As an example, the results of such calculations in the case 
c — ZksIM are presented in Figure 1.1. The dependence of the temperature 
and velocity of the stationary vaporization front on the laser intensity absorbed 
are given in references 1, 14, and 18 for different metals. 

Let us consider now the transient process leading to the formation of a 
stationary evaporation wave. As Equation (1.19) shows, in a stationary wave, 
the spatial scale of temperature distribution near the solid-gas phase bound¬ 
ary is equal to max (x/^4? M -1 )- To derive a very rough estimate of the time 
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required to establish steady-state evaporation at constant absorbed laser in¬ 
tensity labs, we calculate the time of formation of the heated layer near the 
phase boundary. Since in the majority of cases x/K > A* _1 > we see that this 
time is on the order of x/V* 2 - The nature of the process of the establishment 
of steady motion of the phase boundary can be described as follows. 

At the outset the phase boundary is at rest, and the entire laser energy flux 
absorbed is transferred by heat conduction into the bulk of material. Near 
the surface, a temperature gradient of the order I a bs /« is created, which re¬ 
mains roughly constant, as long as the velocity of the boundary is low. Both 
the surface temperature and thickness of the heated layer near the surface 
increase with time as y/t. Due to the very strong temperature dependence of 
the evaporation rate, the evaporation front velocity remains small, in compar¬ 
ison with V 8 , until the temperature attains value very close to T 5 (0). Thus, 
throughout the period that lasts on the order of x/V^ 2 , the phase boundary 
remains practically at rest; for the remainder of the laser-pulse duration, the 
boundary moves with constant velocity V s . It is clear that if the laser-pulse 
length is shorter than x/V^ 2 , then there is no appreciable evaporation of the 
material at all. Hence, we can estimate the critical laser intensity at which 
significant evaporation begins to be: 

labs - pLmVx/n 

where 77 is the laser pulse duration. A similar qualitative analysis shows that 
in the case of x/^s /i -1 , the time required to reach a steady-state evap¬ 
oration regime is on the order of (/iV^) -1 . A more detailed analysis of the 
transient process leading to the formation of the stationary evaporation wave 
is presented in references 14 and 19-21. It is interesting to note that under 
certain conditions, oscillatory changes of the surface temperature and vapor¬ 
ization front velocity during the transient process may be observed. 21 This 
oscillatory behavior is the result of the volume nature of the laser radiation 
absorption, and it disappears if the absorption takes place at the surface. 

We implicitly assumed that the surface reflectivity is constant. The effect 
of surface reflectivity change on the characteristics of the steady-state evap¬ 
oration wave and its formation time was considered by Anisimov et al. 20 We 
shall not specify these calculations 20 because they were derived for the model 
dependence of surface reflectivity on surface temperature. A more realistic 
analysis requires numerical calculations, which are described in Section 1.4. 


1.4 Numerical simulation of laser-produced 
melting and evaporation 

In the previous sections, we considered laser-produced melting and vaporiza¬ 
tion as separate and independent processes. Generally, these processes are 
coupled and should be considered jointly. This circumstance strongly com¬ 
plicates the analysis of the laser-induced phase transitions. Additional com¬ 
plications arise due to laser light reflection from the evaporation front. The 
reflectivity depends on the temperature and phase composition of the surface 
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layer of the material. This layer is strongly nonuniform, and its reflectivity 
can be calculated using Maxwell’s equations. It is clear that a quantitative 
consideration of this problem demands numerical calculations. We shall now 
discuss some of the results of such calculations . 10 

The model employed by Anisimov and Barsukov 10 includes the heat con¬ 
duction equations for the temperatures of solid and liquid phases, with bound¬ 
ary conditions (1.14) and (1.15) on the melting front and (1.18) on the vapor¬ 
ization front. These equations can be written as: 

p s c s dT s /dt = d/dz(K s dT s /dz) + Q 5 (z, t) Z m (t) < z < oo 
picidTi/dt = d/dz(mdTi/dz ) + Qi(z,t) Z v (t) < z < Z m (t) 

The boundary conditions read: 


LmZ m (t) — K s dT s /dz \ z -z m K idTi/dz\z = z m 

L v Z v (t) = K,idTi/dz\ z= z m 23) 

Z v (t) = Vb exp {-U/T v ) T v = T t (Z v ,t) 

Zm(t) = V m c s (Tf -T m )/L m Tf = T s (Z m ,t) 

Here, the subscripts s and / refer to the solid and liquid phases, respectively. 
The laser energy released in the liquid is written as: 

Ql(z,t) = ai(z,t)\E(z,t)\ 2 (1-24) 


where 07 (z, t) is the local electrical conductivity of the liquid material (at laser 
frequency u;o), and E(z,t) is the local amplitude of the electric field in the 
laser beam, its value at the surface being related to the incident laser radiation 
intensity, Jo, by: 

\E(Z v ,t)\ 2 = 47t/ 0 /co 

Here, Co is the speed of light. An equation similar to (1.24) also can be written 
for the solid phase. The electric field inside the metal is described by the wave 
equation: 

d 2 E/dz 2 + kle(z)E = 0 k 0 = w 0 /co (1.25) 

which follows directly from Maxwell’s equations. 

Equations ( 1 . 22 ), (1.23), and (1.25) were solved numerically. The dielectric 
constant of the metal, 6 , was considered dependent on the local temperature 
and phase state of the metal. The temperature dependence of e was cal¬ 
culated using the Drude-Zener model and experimental data on the optical 
properties of metals . 22 Note that the experimentally measured reflectivity of 
some metals may significantly differ from the Drude-Zener reflectivity . 23 This 
difference, however, decreases at high temperatures and low photon energies 
(< 1 eV). The heat conduction equations (1.22) were solved using the front¬ 
capturing numerical scheme, which provides a precise calculation of melting 
and vaporization front positions. Some of the results of these calculations are 
presented in figures 1.2 through 1.6. The results are given for silver; the laser 
intensity equals 10 9 W/cm 2 . The initial reflectivity of cold silver corresponds 
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Figure 1.2: The reflectivity of a silver target as a function of time; 
I = 10 8 W/cm 2 . 


to a high-quality optical surface (R ~ 0.99). In Figure 1.2, the change of sur¬ 
face reflectivity over time is shown. In the absence of melting, the reflectivity 
slowly decreases over time. When the liquid phase appears, the reflectivity 
changes drastically from its “solid” to “liquid” value; but there is no jump in 
reflectivity. Figures 1.3 and 1.4 show the temperature over time at the melting 
and evaporation fronts, respectively. Note that the melting front temperature 
is well above the melting point of silver, which is T m = 1234 K. In Figure 1.5, 
the melting (curve 1) and vaporization (curve 2) front positions are shown as 



Figure 1.3: Temperature on the solid-liquid phase boundary as a function of 
time. 
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t,ns 

Figure 1.4: Temperature on the liquid-vapor phase boundary as a function of 
time. 


functions of time. The melting front propagates in a stationary regime; but 
for the evaporation front, this regime has yet to be reached. Figure 1.6 shows 
the spatial profile of temperature in the vicinity of the vaporization front. As 
was mentioned in Section 1.3, the temperature has a maximum located inside 
the liquid phase, at a distance on the order of fi~ l from the phase boundary. 
Although the difference between the maximum temperature and the temper¬ 
ature at the phase boundary is rather small, it results in a qualitatively new 
phenomenon: the corrugation instability of the phase boundary. 



t,ns 


Figure 1.5: The positions of melting (curve 1) and vaporization (curve 2) 
fronts as functions of time. 
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Figure 1.6: Spatial distribution of temperature in the liquid near the vapor¬ 
ization front. 
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Chapter 2 

INTERACTION 
OF LASER RADIATION 
WITH DIELECTRICS 


2.1 Optical breakdown 

of transparent dielectrics 

In this chapter we will describe the effects that are produced when the laser 
beam is incident on transparent materials. Typical examples of materials un¬ 
der consideration are wide-band-gap dielectrics and semiconductors. When 
the laser intensity is below a certain threshold value, the material is transpar¬ 
ent and the laser beam passes through the material with no apparent effect. 
But when the intensity is above the threshold, absorption abruptly increases, 
and effects such as melting, evaporation, and shock wave generation can oc¬ 
cur. A violent increase in optical absorption of initially transparent material 
caused by a laser pulse is known as laser-induced breakdown. The laser break¬ 
down is, in fact, a special kind of instability of a transparent dielectric in a 
high-intensity electromagnetic field. 

Many different mechanisms of laser breakdown have been described in the 
literature (see references 1-3). The breakdown of real optical materials fre¬ 
quently is associated with the presence of inclusions, impurities, or structural 
imperfections of the material, in which the local absorption of light is much 
higher than the average absorption over the whole sample. In many cases of 
practical interest, the size of such absorbing inclusions is much less than the 
size of the focusing region. 

In pure materials, on the other hand, the above extrinsic mechanisms of 
breakdown can be eliminated to a large extent. This implies a higher break¬ 
down strength. The limiting strength is determined by intrinsic breakdown 
mechanisms, such as multiphoton ionization and electron avalanche. 

Let us first consider the breakdown of transparent dielectrics initiated by 
absorbing inclusions. 4,5 There are two possible mechanisms of development 
of breakdown in this case. In the presence of inclusions of sufficiently large 
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size, the most probable breakdown mechanism is local melting of material 
near the inclusions and formation of cracks as the result of thermal stresses. 
This mechanism of damage takes place, for example, in platinum-containing 
glasses. 6 It is possible, however, that there is another breakdown mechanism 
due to the fact that heating of the transparent medium in the vicinity of an 
absorbing inclusion leads to some additional absorption, which is equivalent to 
an increase in the size of the inclusion. Under certain conditions, this process 
results in the formation of a thermal wave propagating from the absorbing 
inhomogeneity through the volume of the medium. This wave arises as the 
result of thermal instability of the medium near the absorbing inclusion. 

Consider the temperature field near an isolated absorbing inclusion. Let 
us assume for simplicity that the inclusion is spherical and is separated from 
the surrounding dielectric by a sharp boundary. In a quasi-stationary mode, 
the value of energy, which the inclusion acquires from the laser beam, is equal 
to the energy it transfers to the surrounding medium. We denote by I the 
laser intensity near the inclusion. We then can set the heat flux transferred 
from the inclusion to the surrounding dielectric equal to: 

-«VT| r=fl = a(R)I (2.1) 

where k is the thermal conductivity of dielectric, R is the radius of inclusion, 
and a(R) is the ratio of the absorption cross section of the inclusion to the area 
of its surface, 47rR 2 . The factor a(R) depends on the laser wavelength and 
can be calculated from the solution of the problem of light-wave diffraction at 
the inclusion (e.g., reference 7). In all cases of practical interest, a(R) is less 
than unity. 

As mentioned above, the heating of the medium near the inclusion leads to 
some additional light absorption. We assume that a thermodynamic equilib¬ 
rium is established near the inclusion so that the local absorption coefficient 
is governed by a local temperature. The temperature dependence of the ab¬ 
sorption coefficient of media under consideration usually can be described by 
fi{T) = /ioexp(— E/T), where 2E = E g is the forbidden band width. Thus, 
/i(T) is very strong function of temperature. The steady-state temperature 
distribution T s (r) can be found from the solution of the heat conduction equa¬ 
tion: 

kAT s + Ifi{T s ) = 0 (2.2) 

We will show that the steady-state temperature field, T 5 (r), turns out to be 
unstable if the laser intensity, /, exceeds a certain value. To analyze the 
stability of the solution of boundary-value problem given in Equations (2.2) 
and (2.1), we consider the transient heat-conduction problem and, as usual, 
assume that T(r,t) = T s (t) + 0(r,£), Q(r,t) T s (r). For Q(r,t) we then find 
the linear equation: 

cpdO/dt = kAO + V(T 5 ) (2.3) 

where /i'(T 5 ) = dfi(T s )/dT s . We write the solution of Equation (2.3) as: 

0(r,t) = ^2C\tp\(r)exp(-Xt) 

A 


(2.4) 



CHAPTER 2. INTERACTION WITH DIELECTRICS 


27 


where the functions rp\(r) satisfy the Schrodinger equation: 

+ (e + V (r))V> = 0, V ( r ) — Ih'(T s )/k, e = cp\/n (2.5) 

(the index A is omitted from the function rp for brevity). The boundary 
conditions for Equation (2.5) are: 

ip(oc) = 0 and dip/dr\ r=R = 0 

The second condition means that we consider temperature fluctuations in the 
medium while the thermal state of the inclusion remains constant. If among 
the solutions of Equation (2.5) there are those corresponding to negative values 
of A, the steady-state temperature distribution, T 5 , turns out to be unstable. 
By virtue of the previous arguments for the properties of the function /i(T), 
we conclude that the potential V(r) in the Schrodinger equation (2.5) is large 
at r ~ R and decreases rapidly with increasing r. Thus, we are dealing with 
a standard quantum-mechanical problem of a bound state in a deep narrow 
well. Substituting \p = ip/r into (2.5) and integrating the resulting equation 
term by term, we find the instability condition: 


oo 

J fi\T s (r))dr > k/IR (2.6) 

R 

To evaluate the integral in (2.6), we transform the integration variable r to the 
new variable T s and incorporate the fact that /i'(T 5 ) decreases rapidly with 
T s decreasing. Then, the integration can be performed easily and yields: 

IRfi[T s (R))>K\r s (R)\ 

Using the boundary condition (2.1), we arrive at the instability condition: 

Rfi[T s (R)) > a (2.7) 

The left side of Equation (2.7) is, in order of magnitude, the fraction of inci¬ 
dent laser intensity, which is absorbed by heated dielectric material near the 
inclusion. When this fraction becomes larger than the fraction of incident 
laser intensity absorbed by the inclusion itself, the medium turns out to be 
unstable. In order to obtain the threshold intensity, one must solve the bound¬ 
ary value problem (2.1) and (2.2). The calculation carried out in reference 4 
yields the critical laser intensity 

I* « KE/Ra(R) (2.8) 

Note that in the previous analysis, we assumed that the thermal conductivity 
of material does not depend on temperature. Although the breakdown crite¬ 
rion (2.7) is not associated with this assumption and is of a general nature, 
the value of the critical intensity (2.8) may be affected by the assumption 
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made. The temperature dependence of thermal conductivity of the mate- 
rials under consideration is due mainly to the electron component of heat 
conduction, which becomes the dominant energy transfer mechanism at high 
temperatures. The electron heat conduction plays a dual role in the break¬ 
down process. On the one hand, it lowers the temperature, as well as the 
electron density near the surface of an inclusion; on the other hand, it in¬ 
creases the size of the absorbing halo around an inclusion and, consequently, 
enhances the absorbed energy. Therefore, estimating the effect of the electron 
heat conduction on the threshold of forming an absorbing wave is not a trivial 
task. The influence of the electron thermal conductivity in thermally ionized 
dielectric material on the characteristics of its optical breakdown initiated by 
an absorbing inclusion has been studied in reference 8. It was shown that 
allowing for electron thermal conductivity changes the critical intensity (2.8) 
relatively little. 

As can be seen from Equation (2.8), the threshold of breakdown initi¬ 
ated by an individual inclusion depends substantially on its size (absorption 
cross-section). If the medium contains randomly distributed inclusions of dif¬ 
ferent sizes, the breakdown threshold is no longer a completely defined quan¬ 
tity, and we can speak only on the probability of a breakdown under given 
conditions. Thus, the breakdown of a medium containing absorbing inhomo¬ 
geneities should be described by statistical methods. 

Let the distribution of laser intensity over the focal volume be J(r) = 
Jo<£(r). To each intensity value there corresponds a critical size of inclusion 
R*( r), which is related to J(r) by Equation (2.8). Breakdown inside the 
volume V will occur if the volume contains at least one absorbing inclusion 
whose size exceeds the critical value J2*(r) at the corresponding point. As¬ 
suming that the distribution of inclusions is spatially uniform, it is easy to 
write an expression for the breakdown probability in the volume V: 


^[^( r )] = [1 - exp(-nV)] 1 


1 — exp 



(2.9) 


Here n is the average density of absorbing inclusions; f(R) is the distribution 
function of inclusions in size, normalized to unity; and R*(r) is the function 
related to I( r) by Equation (2.8). In the experimental study of breakdown, 
the spatial profile of the intensity usually remains fixed, and the maximum in¬ 
tensity Jo is varied. The breakdown probability in this case is a function of Jo, 
and using Equation (2.9), one can determine the average value and dispersion 
of the random value of the breakdown threshold. The calculations conducted 
in reference 5 show that, when the number of inclusions in the focal volume is 
large, N = nV 1, the average value of breakdown intensity turns out to be 
strongly dependent on the focal volume V and the asymptotic behavior of the 
distribution function of inclusions f(R) at large R. Assuming that as R —► oo 
the distribution function f(R) falls out as f(R) oc R~ k , the average break¬ 
down intensity, < Jo >, was found to be proportional to V~ l ^ k ~ l \ The last 
dependence is in satisfactory agreement with the experimental data 9 if one 
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assumes the exponent k = 2.5. This “size effect” disappears when the focal 
volume and/or average concentration of inclusions is small, N = nV 1. In 
this case, however, the breakdown usually is not due to absorbing inclusions 
in the medium, but to such intrinsic mechanisms as multiphoton ionization 
and electron avalanche. 

The previous analysis of the instability of transparent dielectric materials 
containing absorbing inclusions was based on study of small temperature per¬ 
turbations. This kind of linear analysis can yield only the instability criterion. 
To study the further development of instability, one should solve a much more 
complicated nonlinear problem. As indicated earlier, the consequence of the 
instability is the occurrence of an absorption wave. The propagation of this 
wave from the inclusion to the volume of the medium leads to the expansion 
of the ionized region inside the insulator that results finally in the absorption 
of a considerable portion of the laser energy and the appearance of macro¬ 
scopic damage. In reference 8 numerical modeling has been employed to find 
the qualitative features of this phenomenon. It was assumed that the laser 
radiation impinges on a glass target containing a small metallic inclusion. The 
temperature field in the inclusion and glass matrix was described by the set 
of heat-conduction equations: 

CiPi&Ti/dt = r~ 2 d/dr(r 2 KidTi/dr) + Qi (2.10) 


where i = 1 for the inclusion and i = 2 for the glass matrix. The set of 
equations (2.10) was solved numerically subject to the following initial and 
boundary conditions: 


Ti(r,0) = T 2 (r,0) = T 2 (oo,t) = T 0 

3Ti/0r|r =o = 0 (Ki&Ti/dr - n 2 dT 2 /dr)\ r=zR = 0 


where R is the radius of inclusion. A small inclusion with a radius much 
smaller than the wavelength of laser radiation was considered. If we ignore the 
“shadow” thrown by the inclusion and suppose that the absorption in the halo 
is weak, we can assume that the temperature field has a spherical symmetry, 
as allowed by Equation (2.10). The dielectric in the vicinity of inclusion 
was assumed to be a partially ionized dense plasma. Transport and optical 
properties of this medium were expressed in terms of the effective collision 
frequency. The details of calculations of Qi and ac* are described in reference 8. 
We shall discuss here some of the results of the modeling. Figure 2.1 shows 
the spatial distribution of the temperature for different times. When the 
intensity is below the breakdown threshold defined by (2.8), the stationary 
temperature profile is established in the vicinity of inclusion (dashed curve). 
However, when the intensity exceeds the breakdown threshold, the absorption 
wave is formed, which propagates from the surface of the inclusion to the 
volume of dielectric (full curves). The size of the ionized halo as a function 
of time is shown in Figure 2.2. For laser intensities below the threshold, the 
size of the heated region saturates after some transient process (curve 1). 
For supercritical laser intensities, the front of thermal wave asymptotically 
approaches the regime with constant propagation velocity (curve 2). The 
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Figure 2.1: Temperature profiles at different times. 

Dashed curve—steady-state temperature profile for I = 2 x 10 6 W/cm 2 < 7*, 
Solid curves—temperature in the thermal wave for I = 2 x 10 7 W/cm 2 > 7*, 
I* = 10 7 W/cm 2 . 

glass matrix parameters used in calculations 8 were chosen to correspond to 
the experimental conditions. 10 The velocity of a steady-state absorption wave 
and the maximum temperature reached in this wave are in agreement with 
the experimental results. 10 

We considered the thermal ionization mechanism, in which absorption re¬ 
sults from the appearance of free carriers in the conduction band of a solid. 
Other mechanisms of additional absorption may play roles in laser-induced 
breakdown of transparent dielectrics. There are a number of thermochemical 
processes resulting from the laser heating of a solid. They include the decom¬ 
position of complex compounds and formation of oxides. If the absorption co¬ 
efficient of reaction products at the laser frequency is significantly high, these 
processes lead to thermochemical instability. This instability plays, for exam¬ 
ple, an important role in laser breakdown of polymers. 11 It was established 



t, *iS 


Figure 2.2: Position of the thermal wavefront as a function of time. 

1—7 = 2 x 10 6 W/cm 2 < /*, 2—7 = 2 x 10 7 W/cm 2 > 7*, I* = 10 7 W/cm 2 . 
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that the thermochemical mechanism involving pyrolysis and the formation of 
highly absorbing soot particles is responsible for the laser-induced damage in 
polymethylmethacrylate under the action of continuous laser radiation . 12,13 
The qualitative picture of instability development and absorption halo for¬ 
mation is similar to that of the thermal ionization mechanism. However, the 
pyrolysis reaction requires a long time 13 and does not play a significant role 
in the case of short laser pulses. 

Now, let us briefly consider the intrinsic mechanisms of breakdown of 
transparent dielectrics: multiphoton and avalanche ionization. The electron 
avalanche is the main mechanism of breakdown of gases, as well as solid 
and liquid dielectrics in both d.c. and high-frequency electric fields (see refer¬ 
ences 14-16). The avalanche is the most probable mechanism of laser break¬ 
down in ultrapure transparent dielectrics, if the laser pulse length is not too 
short and the focal volume is not too small. As in the case of thermal ion¬ 
ization, the electrons in the conduction band acquire energy from the laser 
beam due to electron-phonon collisions. High-energy electrons then can pro¬ 
duce the ionization, leading to electron multiplication. The major difference 
between the thermal and avalanche ionization is in the distribution function 
of electrons in the conduction band. In the case of avalanche, it is strongly 
nonequilibrium and should be found from the solution of the quantum kinetic 
equation . 16 As applied to laser breakdown of transparent dielectrics, different 
approximate solutions of this equation have been obtained in references 17-19. 
The distribution function of electrons usually is assumed to have the form: 

/(P, t) = /(p) exp( 7 t) 

where p is the momentum of the electron, and the avalanche development 
constant 7 is calculated as a function of crystal parameters and laser pulse 
characteristics. The electron distribution function is supposed to be indepen¬ 
dent of spatial coordinates. When 7 is determined, the threshold intensity 
is calculated from the breakdown criterion, which can be taken in the form: 
777 > C, where C is either a constant or a weak (logarithmic) function of 
the crystal parameters and concentration of seed electrons. It generally is 
accepted that the concentration of seed electrons, no, is not too small, so that 
txqV > 1 , where V is the focal volume. In the opposite case, n$V < 1 , the 
breakdown threshold becomes a random quantity determined by the proba¬ 
bility of appearance of a seed electron inside the focal volume. 

Multiphoton ionization plays a dual role in the laser-induced breakdown 
of ultrapure insulators. On one hand, it is an auxiliary process responsible for 
supplying seed electrons necessary for the avalanche breakdown. On the other 
hand, the multiphoton ionization is a competitive mechanism of breakdown, 
which can provide, under certain conditions, a higher rate of generation of 
electrons in the conduction band than in the electron avalanche. The analysis 
of the relative contribution of the two processes was performed by Gorshkov 
et al . 20 It was shown that the dependence of the breakdown thresholds on 
the laser pulse length is different for avalanche and multiphoton breakdown 
mechanisms. Multiphoton ionization proved to be the dominant breakdown 
mechanism in picosecond and subpicosecond ranges of pulse lengths. 
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2.2 Stationary evaporation waves 
in transparent dielectrics 

In Section 2.1 we described the propagation of a laser-supported absorption 
wave inside a transparent dielectric. The formation of the wave was initiated 
by an absorbing inclusion. In this section we shall consider the wave initiated 
by a “seed” located at the surface of transparent material. The experiments 21 
show that in this case the propagation of absorption wave from the seed to the 
volume of dielectric is accompanied by the surface evaporation. Under cer¬ 
tain conditions, synchronous propagation of the absorption and vaporization 
waves appeared to be possible. It also was shown that, for millisecond laser 
pulses, the stationary evaporation regime of the entrance and exit surfaces is 
established in silicate glasses and fused quartz. 21 

Consider a planar, steady-state evaporation wave propagating through the 
material with a strongly temperature-dependent absorption coefficient. Let 
the velocity of the wave be V a . The temperature distribution in the wave is 
described by the heat conduction equation transformed into the moving frame, 
which is bound to the vaporization front: 

cpV s dT s /dz + d/dz[K(T s )dT s /dz] + Q(T S ) = 0 (2.12) 

where Q(T) = I(z)p(T). The equation for laser radiation intensity can be 
written as: 

dl/dz = T#i(T a )J (2.13) 

where minus and plus signs in the right-hand side correspond to the evap¬ 
oration of the entrance and exit surfaces, respectively (see the schematic in 
Figure 2.3). The boundary conditions for Equations (2.12) and (2.13) are: 

K(T s )dT s /dz\ z = 0 = pV s AH « pV s L v T s ( oo) = 0 (2.14) 

1(0) = Jo (evaporation of the entrance surface) 

or J(oo) = Jo (evaporation of the exit surface) 

Note that only a part of the incident laser intensity is absorbed by the di¬ 
electric. The transmitted laser intensity (J(oo) in the case of entrance surface 
evaporation, and 1(0) in the case of exit surface evaporation), is not known 
a priori and should be determined in the course of the solution. As in Sec¬ 
tion 2.1, we assume the dependence p(T) in the form p(T) = /xoexp(— E/T) 
where E = E g / 2. The dependence of vaporization front velocity on the tem¬ 
perature is taken in the form of Equation (1.16): V s = Vbexp[— U/T s (0)\. 

Substituting Equation (2.13) into Equation (2.12) we obtain the first in¬ 
tegral, a form of the energy-flux conservation: 

K(T s )dT s /dz = ±[J - J(oo)] - cpV s T s (2.15) 

If we set z = 0 in Equation (2.15) and employ the boundary conditions (2.14), 
we obtain the familiar formula: 


P V s (L v + CT S ( 0)) = ±[J(0) - J(OO)] = labs 


( 2 . 16 ) 
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Figure 2.3: Schematic of the experiment for observing the evaporation of 
entrance (1) and exit (2) surfaces of a transparent dielectric. Qualitative 
profiles of the temperature and laser intensity are shown for both cases. 


which coincides, formally, with (1.20) and can be combined with (1.16) to 
calculate the velocity V s and surface temperature T 3 (0). In practice, however, 
this cannot be done since the absorbed laser intensity I a bs is not known. 

The determination of temperature profile in the stationary evaporation 
wave is reduced to the integration of the coupled set of equations (2.13) and 
(2.15), subject to the boundary conditions 

T s { oo)=0, 1(0) = I 0 (2.17) 

(or I( oo) = Iq for exit surface evaporation). The velocity V 3 is an eigenvalue 
of problem (2.13) and (2.15). It can be shown that, due to the invariance of 
the above problem under the transformation z —> z 4* zo, z 0 = constant, the 
boundary conditions (2.17) are sufficient to determine the solution T s (z) and 
eigenvalue V s uniquely. 

Equations (2.13) and (2.15) were solved numerically for fused quartz by 
Aleshin et al. 21 Figure 2.4 shows the dependence of the fraction of laser energy 
absorbed in the target, A = I a bs/lo, on the incident light intensity Jo- Curves 1 
and 2 correspond to the evaporation of incident and exit surface, respectively. 
The absorption decreases as the laser intensity increases in the range 10 6 - 
10 7 W/cm 2 and remains approximately constant for intensities in the range 
10 7 -10 8 W/cm 2 . The absorption is higher for the evaporation of exit surface 
than for the evaporation of entrance surface. At high intensities, when the 
absorbing layer becomes optically thin, the spatial profile of laser intensity 
is independent of the direction in which the radiation propagates, and the 
difference in the evaporation of the entrance and exit surfaces disappears. 
The stationary velocity of evaporation wave V s is shown in Figure 2.5 as 
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Figure 2.4: Absorption A = I a bs/I o of quartz glass as a function of incident 
laser intensity: 1—entrance surface evaporation, 2—exit surface evaporation. 


a function of the incident light intensity. The velocity is proportional to the 
laser intensity in the range Iq > 10 7 W/cm 2 , where the absorption is constant. 

The calculations show that in the case of exit surface evaporation, a sta¬ 
tionary solution exists if the absorption is lower than a certain critical value 
A < A* (or the incident intensity is higher than some threshold value I *). 
When A > A*, the velocity of the absorption wave becomes higher than the 
velocity of the evaporation front. This results in the separation of the absorp¬ 
tion layer from the phase boundary and in the dumping of vaporization. Thus 
we see that, under certain conditions, the laser evaporation of the exit surface 



Figure 2.5: Stationary velocity of the evaporation wave as a function of inci¬ 
dent laser intensity: 1—entrance surface evaporation, 2—exit surface evapo¬ 
ration. 
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of a transparent dielectric can be unstable. This instability will be studied in 
Chapter 6. In the case of entrance surface evaporation, a stationary solution 
exists for any intensity value of the incident radiation. 

Since the absorption coefficient /i(T) is a strong function of temperature, 
the heat source in Equation (2.12) is localized in a thin layer near the max¬ 
imum of temperature. This fact has been used for approximate solution of 
Equations (2.12) and (2.13). 22 The results obtained by Kovalev 22 are in a 
good agreement with the numerical solution 21 at fairly high laser intensities 
when the absorption layer is optically thin. 


2.3 UV-laser ablation of polymers 

Material removal caused by laser irradiation is often denoted as laser ablation. 
Recent experiments 23,24 show that for polymer materials the kinetics and dy¬ 
namics of material removal strongly differ for UV and IR laser radiation. This 
means that the mechanism of laser-polymer interaction depends on the photon 
energy, hv. Thus, the traditional question of the competition of thermal and 
nonthermal mechanisms in laser ablation arises. A similar question also has 
been discussed in the laser annealing of semiconductors. 25 Recently, many ex¬ 
periments on UV-laser ablation of polymers have been performed with excimer 
lasers and pulse durations on the order of 10-30 ns. These experiments show 
that the products of ablation of organic polymers consist of low-molecular 
fragments (carbon, simple organic molecules) and medium-molecular frag¬ 
ments (monomers and parts of them). 26 To create such fragments, covalent 
bonds within the original polymer must be broken. If these bonds are broken 
due to thermal motion of atoms in molecules, the ablation velocity can be 
described as: 

V = V 0 exp(-E/T 0 ) (2.18) 

where Tq is the surface temperature and E is the bond-breaking energy. Typ¬ 
ically, this energy is about 3-5 eV. Thus, experimentally observed ablation 
velocities of the order of 0.1 /xm/pulse « 1 m/s can be obtained only at 
very high temperatures. The theoretical model 27 yields temperatures about 
To ~(6-12)xl0 3 K near the ablation threshold. These values contradict the 
experimental data on the vibrational temperatures of ablation products, 26 
T V ib ~(2-3)xl0 3 K. The experiments also show that the UV-laser ablation 
of heat-sensitive polymers can be performed without damaging the remaining 
material, specifically with no trace of melting. The above facts cannot be 
explained on the basis of the thermal mechanism of ablation. 

Photochemical ablation deals with direct breaking of chemical bonds as the 
result of light absorption. For the case of one-photon absorption this yields 
the dependence of the ablated layer thickness A/i on the laser fluence $ in the 
form: 

Ah = /x -1 log($/$th) (2.19) 

where $>th is the threshold fluence for ablation. In recent experiments, 28 the 
dependence h($) has been measured for polyimide with four different laser 
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wavelengths. Only with the ArF-laser (A = 193 nm) was agreement obtained 
with Equation (2.19) at $ « All other wavelengths revealed dependen¬ 
cies that disagree with (2.19). Thus, the experimental results can neither be 
explained on the basis of a purely photochemical process, nor on the basis of 
a purely thermal process. They can, however, be described on the basis of a 
photophysical ablation mechanism. 29,30 The main idea of this mechanism is 
that the removal of electronically excited species from the surface requires a 
smaller activation energy than the removal of species in the electronic ground 
state. 

Let us assume that the total ablation velocity is determined by both 
ground-state species A and excited-state species A*. The excitation process 
A A* is induced by one-photon absorption. The dissipation of the excita¬ 
tion energy A* —> A occurs either via stimulated emission, or via nonradiative 
transitions. The nonradiative transitions will be characterized by a single 
thermal relaxation time r^, which is assumed to be independent of tempera¬ 
ture. We assume also that a part of the excitation energy is lost via activated 
desorption of excited species. 

Consider a planar ablation front. Neglecting any cooperative phenomena 
in the desorption process, the ablation front velocity can be described as: 


V = xVf exp(—E*/To) 4- (1 - x)V 0 exp(-E/T 0 ) 

where x = N* /(N* 4- N) is the normalized density of excited species; N and 
N * are the number densities of ground-state and excited-state species, respec¬ 
tively; and E and E* are their respective activation energies of desorption. 
We assume E* E and, additionally, that E is on the order on the bond¬ 
breaking energy. In a coordinate frame connected with the ablation front, the 
density of species can be described by: 

ON*/dt = VdN*/dz + {aI/hv){N - N*) - N*/t t 
dN/dt = VdN/dz - (oI/hv)(N - N*) + N*/t t 

where a is the absorption cross section. Note that N 4- N* = iVo = constant 
is the total number density of chromophores. According to Sauerbrey and 
Pettit, 31 it is typically Nq ~ 6 x 10 21 cm -3 . 

The propagation of the laser light is described by the equation: 

dl/dz = -a(N - N*)I (2.21) 

The energy transfer in the solid polymer is described by the heat-conduction 
equation: 

8T/dt = VdT/dz + x&T/dz 2 + (x/»e)Q (2.22) 

where n is the thermal conductivity, which is assumed to be temperature- 
independent, and x is the thermal diffusivity. The heat source Q is due to the 
nonradiative transitions and is given by: 


( 2 . 20 ) 


Q — hvN* /tt 


( 2 . 23 ) 
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The boundary conditions at the ablation front z = 0 can be written as 
KdT/dz = p[xAH*V* exp(—E*/T) -f (1 - x)AHV 0 exp (-E/T)} 

(2.24) 

7(0, t) — 7o at 2 — 0 

In addition to (2.24), we employ the boundary conditions 

iV(z, t) —> iV 0 , N*(z, t) —> 0, T(z, £) —> 0 as z —► oo (2.25) 

We now concentrate on stationary solutions of the boundary value problem 
(2.20), (2.21), (2.22), (2.24), and (2.25). If we set all time derivatives equal 
to zero in Equations (2.20)-(2.22) and assume that 7o = constant, we obtain 
the first integral 

I s = V s (cpT s + hvN* s ) 4- ndTs/dz (2.26) 

which represents the energy-flux conservation. Here, the subscript s indicates 
stationary quantities. Combining Equation (2.26) at z = 0 with boundary 
condition (2.24), we obtain the equation connecting the number density of 
excited-state species N*( 0) with the surface temperature To = T s ( 0); it can 
be written as: 

Axl -f Bx o 4* C = 0 (2.27) 

where A = hvN 0 [V 0 exp(-£/T 0 ) - VJ exp(-£7T 0 )] 

B = (cpT 0 /hvN 0 )A 4- p [AHV 0 exp{-E/T 0 ) - AH*VJ exp{-E*/T 0 )} 

C = 7 0 - (N 0 hv 4- cpT 0 4- pAH)V 0 exp (~E/T 0 ) 
and x 0 = N*(0 )/Nq 

With tt = constant, Equations (2.20) and (2.21) become independent 
from (2.22), and reduce to the single first-order differential equation 29 

V s dN*/dl s = 1 jhv - N*/t T (tI s {No - 2JV 5 *) (2.28) 

with the boundary condition: N* = 0 at 7 = 0 (i.e., when z —► oo). Integration 
of (2.28) from I s = 0 to I s = 7o permits the calculation of N*( 0) = xoNo. We 
then have the set of equations (2.27) and (2.28) to calculate the parameters of 
a stationary ablation wave: N*( 0), To, and V a . Numerical calculation of these 
parameters was performed in reference 29. Figure 2.6 shows the normalized 
ablation velocity V = VshvNo/Ib as a function of normalized laser intensity, 
Io/h for different relaxation times tt- The value of constant h was chosen 
to be 10 7 W/cm 2 . As seen from Figure 2.6, the ablation velocity increases as 
the relaxation time increases. 

The dependence of the surface temperature on the laser intensity is plotted 
in Figure 2.7 for different values of tt . The full curves refer to the activation 
energy E = 3 eV and E* = 0.3 eV. The dashed curves have been calculated 
for E = 4.5 eV. The temperature is given in the units T*, = 4 x 10 3 K. We 
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Figure 2.6: Normalized ablation velocity V = V s Nohv/Ib as a function of laser 
intensity at different relaxation times vp. 

I*, = 10 7 W/cm 2 , Tp in units 10~ 8 s. 



lo/lb 


Figure 2.7: Ablation front temperature as a function of laser intensity for 
different rp. 

lb = 10 7 W/cm 2 , Tb — 4 x 10 3 K, rp —in units 10" 8 s. 
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Figure 2.8: Concentration of excited species as a function of laser intensity. 
I b = 10 7 , rp —in units 10” 8 s. 


see that the surface temperature decreases as the result of energy relaxation. 
The “cold” ablation related to the deposition of excited species becomes in¬ 
creasingly important. With higher activation energies, the difference between 
thermal and photophysical ablation becomes more pronounced. In Figure 2.8 
the concentration of excited species N*( 0) is given as a function of laser in¬ 
tensity at different relaxation times. Note that N* saturates at high laser 
intensities. 

The spatial distribution of the normalized temperature T 5 (z)/To, is shown 
in Figure 2.9 for different values of Jo and rp. The coordinate z is given 
in units zo = (crNo)" 1 . The width of the temperature distribution increases 
with intensity Jo and relaxation time rp. On an expanded spatial scale, a 
maximum of temperature at some finite distance z C zo can be observed. 
This maximum is related to the finite penetration depth of the incident laser 
radiation. The spatial profile of N*(z)/Nq is shown in Figure 2.10 for different 
Jo and rp. The concentration of excited species reaches its maximum at the 
surface z = 0. 

The model described above explains many features observed in UV-laser 
ablation of organic polymers. With activation energies of E = 3 eV to 5 
eV, E * = 0.3 eV and relaxation times rp — 1 ns the model predicts real¬ 
istic ablation rates at moderate surface temperatures of about 2000 K. The 
model also predicts that, although the stationary temperature has a max¬ 
imum at a finite depth, the planar ablation wave is stable with respect to 
small perturbations. 30 This explains why pulsed UV-laser ablation provides 
smooth surfaces with polymers and rough surfaces with metals. 
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Figure 2.9: Spatial temperature profiles at different Jo and vr- 



Figure 2.10: Spatial profiles of the concentration of excited species for different 
Jo and tt • 
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Chapter 3 


PROCESSES 

IN THE VAPOR PLUME 


In this chapter we shall consider the processes that take place in the laser- 
produced vapor cloud. We shall study the structure of a vapor flow which is 
important for the analysis of the interaction of incident laser radiation with the 
evaporated material. We also shall discuss the effect of the ambient gaseous 
atmosphere on vapor-cloud dynamics, including the instability of the contact 
boundary between the expanding vapor and ambient gas. 


3.1 Hydrodynamic boundary conditions 
for strong evaporation 

Consider the steady-state evaporation of a solid (liquid) into a vacuum. Let 
the temperature of the surface of the condensed phase be T s . Let us choose a 
coordinate frame connected with the vaporization front with z-axis directed 
from the condensed phase to the gaseous one. It generally is accepted (ref¬ 
erences 1 and 2) that the atoms emitted from the surface of the condensed 
phase have the Maxwellian velocity distribution for v z > 0, with number den¬ 
sity equal to the equilibrium number density of the saturated vapor n s (T s ). 
For v z < 0, a good approximation for the distribution function is / « 0. 
We therefore see that the velocity distribution of atoms near the surface is 
substantially different from the local equilibrium distribution. Consequently, 
in the immediate vicinity of the vaporizing surface, there is a layer of sev¬ 
eral mean free paths length, in which the distribution approaches equilibrium. 
This so-called Knudsen layer must be considered a discontinuity in the hy¬ 
drodynamic treatment. The kinetic equation must be solved to determine the 
structure of this region and the values of hydrodynamic variables beyond the 
discontinuity. This analysis was conducted by Anisimov 3 . 

The problem under consideration resembles that of a strong shock wave (or 
rather a “strong rarefaction wave”): just as in a strong shock, a considerable 
change of the distribution function occurs, on the scale of several mean free 
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paths. We can use, therefore, the approach proposed by Tamm 4 and Mott- 
Smith 5 in their studies of strong shock-wave structure. The special feature 
of this method is an approximation of the distribution function within the 
discontinuity region by the sum of distribution functions before and after 
the discontinuity with coordinate-dependent coefficients. This approach is 
acceptable if the region of the main change of distribution function is narrow. 

The distribution function in the Knudsen layer can be written in the form: 


where 

/(z, v) = a(z)/i(v) + [1 - a(z)]/ 2 (v) 

(3.1) 


r n s (T 3 )(M/2nkT s ) 3 / 2 exp{-Mv 2 /2kT s ) 

c* 

N 

V 

0 

/l(v) = 

| 0/2 (v) 

c* 

N 

A 

0 


/ 2 (v) = ni(M/2nkTi) 3/2 exp{-(M/2kTi)[(v z - ui) 2 + vl + *$} 


Here a(z) is the unknown function that satisfies the conditions a( 0 ) = 1 
and a(oo) = 0 , T s is the surface temperature, and n s (T s ) is the saturated 
vapor density at this temperature. The values of ni, ui, and T\ refer to an 
equilibrium state beyond the discontinuity and are related in the steady-state 
evaporation regime by the Jouguet condition : 3 


ui=c s {T{) (3.2) 

where c s (T) = ( 7 fc^T/M ) 1 / 2 is the speed of sound. The last condition follows 
from the fact that the expansion of vapor into a vacuum is described by the 
centered rarefaction wave . 6 

The conservation laws of mass, momentum, and energy fluxes hold within 
the discontinuity region: 


J dvv z f(z,v) = Ci 

J dvv 2 z f{z,v) = C 2 (3.3) 

J dvv z v 2 f(z,V) = Cz 

Equations (3.2) and (3.3) form the set of coupled equations for the calcula¬ 
tion of flow parameters at the sonic point and the constant (3. Solving these 
equation, we obtain : 3 

/3 = 6.29, T\ = 0.67T 5 , m = 0.31n 5 (T 5 ) (3.4) 

Note that for the adiabatic index 7 , the value 5/3 was assumed, corresponding 
to a monatomic gas, since the equilibrium of internal degrees of freedom is 
established at a distance from the phase boundary that is much larger than 
the mean free path of atoms. 

It can be shown using Equations (3.4) that the flux j _ of atoms condensed 
back on the phase boundary amounts to approximately 0.18 of the flux of evap¬ 
orated atoms = (fcBT 5 /27rM) 1 / 2 n 5 (T 5 ), i.e., We now can obtain 
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more precise formulas than those derived in Chapter 1 for the vaporization 
front velocity, vapor temperature, condensed phase temperature, and other 
parameters. Below are some results presented without detailed calculations. 
The surface temperature T s = T s { 0) in the stationary evaporation wave is 
determined by the equation (see Equations (1.20) and (1.21) for comparison): 

n a (T s )(Mk B T s ) 1/2 (L v + 2.2 k B T s /M) = 3.1 I aba 


where L v is the latent heat of vaporization per unit mass and I a b S is the laser 
intensity absorbed. The velocity of the vaporization front is determined by 
the formula: 



labs 

p(L v 4- 2.2fc#T s /M) 


These results differ slightly from those derived without considering the vapor 
expansion dynamics. It is important to note that under certain conditions, 
the vapor at the sonic point appears to be supersaturated. This means that 
the condensation may begin very near the evaporating surface. To estimate 
the condition for the supersaturation to occur, we compare the vapor den¬ 
sity at the sonic point n\ = 0.31n s (T s ) with the saturated vapor density 
corresponding to the temperature T\ = 0.67T S . Assuming the exponential 
temperature dependence of saturated vapor density, we obtain the condition: 
T s < 0AML v /kB- We thus see that at moderate laser intensities, when the 
surface temperature is considerably lower than U = ML v /ks, the vapor near 
the surface is expected to be in the supersaturated state. 

The foregoing considered the expansion of vapor into vacuum, but the same 
method of derivation of hydrodynamic boundary conditions can be employed 
in the case of expansion of the vapor into an ambient gaseous atmosphere. 
In this case the hydrodynamic velocity of vapor at the outer boundary of the 
Knudsen layer is determined by the external flow and is no longer equal to 
the local sound velocity. The jumps of density and temperature across the 
Knudsen layer become functions of the Mach number M\ of the external flow 
at the outer boundary of the Knudsen layer. To determine the boundary 
conditions, we have to modify Equation (3.2) as follows: 


u\ = Mic s (Ti) (3.5) 

and solve the set of equations (3.3) and (3.5). The dependences Ti(M\) and 
n i(Mi) resulting from the numerical solution of (3.3) and (3.5), can be ap¬ 
proximated by the following analytical formulas: 


Ti «r 4 (l -0.33Mi), n x « n s (T s )/(l + 2.2Mj) 

Again, we can compare the vapor density n i on the outer boundary of the 
Knudsen layer with the saturated vapor density corresponding to the tempera¬ 
ture T \. The comparison shows that the condition required for supersaturation 
to occur depends on the Mach number M\. For example, at M\ = 0.5 we have 
the condition T s < 0.2 L v M/ks, and at M\ = 0.25 T s < 0.1 6L v M/Hb- Thus, 
the range of surface temperatures where supersaturation occurs decreases as 
Mi decreases. 
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The above hydrodynamic boundary conditions are similar (if we extend the 
analogy with the problem of the strong shock wave) to the Hugoniot relations. 
The structure of vapor flow in the Knudsen layer (similar to the internal 
structure of a shock wave) can be found from the solution of the Boltzmann 
kinetic equation. Numerical methods have been employed widely to solve 
this equation for both stationary (references 7-9 and references therein) and 
nonstationary 10 evaporation. We will not discuss these works here. Notice 
only that the hydrodynamic boundary conditions resulting from the numerical 
solution of the Boltzmann equation are in close agreement with the previous 
simple analysis. 


3.2 Condensation in expanding vapor 

It was shown in the previous section that, under certain conditions, the vapor 
on the outer boundary of the Knudsen layer is in a supersaturated state. This 
state is unstable, so the condensation of the vapor will take place immediately 
beyond the Knudsen layer. 3,11 In typical experimental conditions ( T s U), 
the supersaturation is high; the condensation therefore occurs very rapidly 
in a very narrow region, which can be regarded as a surface of discontinuity 
(condensation shock) separating the original vapor from a vapor containing 
droplets of condensed phase. 6 The vapor parameters before and after the dis¬ 
continuity are related by the laws of conservation of mass, momentum, and 
energy fluxes. For one-dimensional flow, they read: 

P\U\ = P 2 U 2 

Pi + P\u\ = P 2 4 - P 2 U 2 ( 3 . 6 ) 

Hi -f uJ/2 = H 2 4" u 2 /2 

Here subscripts 1 and 2 refer to the states on the outer boundary of the 
Knudsen layer and on the outer boundary of condensation shock, respectively. 
The density of equilibrium two-phase system is p = Mn, n = n c -f n s , where 
n s is the saturated vapor number density and n c is the number of atoms in 
condensed phase per unit volume of the two-phase system. Introducing the 
degree of condensation, a = n c /n, we can write the pressure and internal 
energy of the two-phase system in the following form: 

p = p s (T) = n( 1 - a)k,BT, e = (1 - a)L v + 3(1 + a)ksT/2M ( 3 . 7 ) 

Here we assumed that the heat capacity of the condensed phase is 3k b per 
atom. Using the condition for the adiabaticity dS = de + pd(l/p) = 0 and the 
expression for the saturated vapor pressure: 2 

P.(T) = B 0 T~ 1 / 2 exp(-[//T) 

we can obtain the equation for the adiabatic of the saturated vapor in the 
form: 


(y - l/2)da + (3 /y -fa - 1 )dy = 0 
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where y = U/T. The solution of this equation, subject to the initial condition 
a = c *2 when y = y 2 , is of the form: 


(1 - a)(y - 1/2) = (1 - a 2 )(y 2 - 1/2) + 3 log (y/y 2 ) (3.8) 


Using Equations (3.7) and (3.8), we obtain, after some straightforward calcu¬ 
lations, a formula for the sound speed in two-phase mixture: 


C, = (dp/dp)l /2 = 


{k B T/M)W(l-a)(y-l/2) 
[(y 2 — 3y + 3/4)(l — a) 4- 3] 1/2 


(3.9) 


Since the above analysis is justified when T <C ML v /fc B , Equation (3.9) can be 

approximately reduced to c s « [(1 — aJ/c^T/M] 1 ^ 2 . Finally, using Equation 
(3.7) we can rewrite the third equation (3.6) in the form: 

5A: b Xi /2A/ -f- L v -f- u?/2 

(3.10) 

= (1 - a 2 )(5fc B T 2 /2M 4- L v ) 4- 3 a 2 k B T 2 /M 4- u 2 2 /2 V ' 


In the presence of condensation shock, the Jouguet condition (3.2) should be 
set on the outer boundary of the condensation shock: 


U 2 = c 5 (T 2 ,a 2 ) (3.11) 

where c s is given by Equation (3.9). 

Let us turn now to Equations (3.6). According to the results of Section 3.1, 
we can express the left sides of these equations as functions of a single param¬ 
eter, u\. Therefore Equations (3.6), supplemented by the Jouguet condition 
in the form (3.11), constitute a set of equations to determine three unknown 
parameters: T 2 , a 2 , and u\. To simplify the calculations, we suppose that 
T s -C MLy/ks • In this case, as one can see from (3.10), a 2 is of the order 
of a s a ksTg/MLy < 1, so we can neglect a 2 in the first two equations 
(3.6) and then use Equation (3.10) to calculate a 2 . After simple calculations, 
we obtain the hydrodynamic boundary conditions on the outer boundary of 
condensation discontinuity, in the following form: 

T 2 /T 8 « 1 n 2 /n s (T s ) « 0.28 

(3.12) 

v>2 « ( ksTs/M ) a 2 « 0.65fc B T 3 /ML v 


3.3 Dynamics of vapor expansion: 
one-dimensional flow 

The motion of vapor beyond the condensation shock obeys the gas-dynamics 
equations. Since the vapor is transparent, its motion is isentropic and is de¬ 
scribed by the self-similar solution of gas-dynamics equations. 6 For the case 
of one-dimensional expansion of a vapor, this solution depends on the sim¬ 
ilarity variable £ = z/t, where the z-axis is directed perpendicularly to the 
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phase boundary. Introducing the similarity variable into the equations of gas 
dynamics, we obtain: 6 


(u - £)dp/d£ + pdu/d£ = 0 
{u-(,)pdu/d(, + c%dp/d£ = 0 


(3.13) 


One must add to the set of equations (3.13) a condition for the adiabaticity 
of an equilibrium two-phase flow. We write this condition in the form (3.8). 
From Equations (3.13) we obtain, after standard transformations: 6,11 

P2 

u = £ + c s = u 2 -f J c s dp/p (3-14) 

p 

Using Equations (3.8) and (3.9), we obtain from Equation (3.14) in the limiting 
case ksTg/M L v : 


u(0/c 2 « l-([//T 2 )log[T(0/T 2 ] 
* £/c 2 +T(0/r 2 


(3.15) 


Equations (3.15) define the velocity and temperature of vapor as functions 
of the similarity variable £. We see that the dependence of the velocity on £ 
is approximately linear, and the temperature decreases exponentially with £ 
increasing. 

It should be noted that the solution (3.15) is valid only in a limited range 
of £. Adiabatic expansion of the two-phase mixture results in a very fast drop 
in vapor density and, hence, in the “freezing” of condensation. 12 At some 
density, the degree of condensation a reaches its maximum value a max and 
remains constant during subsequent expansion. A more detailed analysis of 
the freezing of condensation is described in reference 11. It was shown that, 
under certain conditions, the freezing can occur in the immediate vicinity of 
the phase boundary. The vapor expansion is described in this case by the 
simple ideal-gas relations given in references 6 and 12. 

We have considered the expansion of vapor into vacuum. This considera¬ 
tion remains approximately valid when the pressure of ambient gas is much 
smaller than the saturated vapor pressure near the vaporization front. If, 
however, these two pressures are not very different, the effect of ambient at¬ 
mosphere must be taken into account. The structure of one-dimensional flow 
in this case is shown schematically in Figure 3.1. The expanding vapor succes¬ 
sively passes through a (1) Knudsen layer, in which the velocity distribution 
of evaporated atoms acquires the equilibrium Maxwellian form; (2) a conden¬ 
sation discontinuity, in which the supersaturated vapor is partially condensed; 
(3) a rarefaction wave; and, (4) a uniform flow region. The outer boundary of 
this region is (5) the contact surface, which separates the evaporated material 
and the gaseous ambient. The evaporated material acts on the surrounding 
gas much as a piston generating (7) a shock wave in the gas. The shock- 
compressed gas occupies (6) a uniform flow region, between the shock wave 
and contact surface. The solutions of gas-dynamics equations describing the 



CHAPTER 3 . PROCESSES IN THE VAPOR PLUME 


51 



Figure 3.1: Expansion of vapor into a gaseous atmosphere. Qualitative struc¬ 
ture of flow. 1 —Knudsen layer; 2—condensation discontinuity; 3—rarefaction 
wave; 4—region of uniform vapor flow; 5—contact boundary; 6 —uniform flow 
of shock-compressed gas; 7—shock wave. 


flow characteristics in regions (3), (4), and ( 6 ) are well known . 6,12 These so¬ 
lutions should satisfy the boundary conditions on the discontinuity surfaces 
(2), (5), and (7). The boundary conditions on the shock wave follow from the 
conservation of mass, momentum, and energy fluxes across the shock front. 
These conditions can be written as : 6 

u 6 = 2D(1 - M 0 ~ 2 )/(7 + 1), Pq = Pa [l + 27 (Mq - l )/(7 + 1)] (3.16) 

where D is the shock wave velocity, Mq = D/c a is the Mach number of the 
shock wave, c a is the sound speed in the ambient gas, p a is the external 
pressure and 7 is the adiabatic exponent of the gas. The conditions on the 
contact boundary (5) are: 


uq = U 4 and pe = p± (3.17) 

Combining the boundary conditions discussed above on the condensation 
shock with Equations (3.14), (3.16), and (3.17), one can determine all charac¬ 
teristics of gas and vapor flows as functions of the surface temperature T s and 
external pressure p a . Numerical calculations for aluminum have been carried 
out by Igoshin and Kurochkin . 13 On the basis of these calculations, approx¬ 
imate analytical equations were proposed for making numerical estimates of 
the flow parameters. 


3.4 Dynamics of vapor expansion: 

spatial structure of the vapor plume 

Analysis of vapor expansion dynamics using one-dimensional approximation 
is justified when the thickness of the vapor layer A z oc c s ri is much smaller 
than the laser spot size. This condition usually is satisfied for nanosecond 
laser pulses. The formation of the vapor/plasma cloud near the surface of 
a target can be described in terms of the one-dimensional model. However, 
the subsequent expansion of the cloud into a vacuum or ambient gas is a 
two- or three-dimensional process. The dynamics of this process have been 
studied experimentally for different laser sources and target materials . 14-18 
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In this section we will consider a simple analytical model that yields a correct 
qualitative explanation of the basic experimental results. The model is based 
on a particular solution of the gas dynamics equations that results from the 
invariance of these equations under a Lie group transformation. 19 We suppose 
that the initial vapor cloud formed near the target surface has an ellipsoidal 
form (with semiaxes Xq , Yo> and Zq) and consider its expansion into vacuum. 
The expansion of the cloud is described by the gas dynamics equations: 

dp/dt + div(pu) = 0 

du/dt + (uV)u + (l/p)Vp = 0 (3.18) 

dS/dt -f u VS = 0 

We assume that the vapor can be described by the ideal gas equation of state 
with the adiabatic exponent 7. 

Consider a fluid motion in which the coordinates of a fluid particle r^t), 
(i = x,y,z), change according to the equation: 20 

ri (t)=]TF ifc (t)r fc ( 0) (3.19) 

k 


where r^(0) are the initial coordinates. Transformation (3.19) describes uni¬ 
form deformations and rotations of a fluid mass. If we ignore the rotation, the 
matrix Fik(t) is reduced to a diagonal form: 

X(t)/X 0 0 0 

Fifc(t) = 0 Y(t)/Y 0 0 (3.20) 

0 0 Z(t)/Z 0 

Equation (3.19) with the matrix given by (3.20) describes the conversion 
of fluid ellipsoid with semiaxes Xo, Vo? and Zq into another ellipsoid, with 
semiaxes X, T, and Z. According to Equation (3.19), the fluid velocity at 
the point (rr, y, z) is proportional to the radius vector of this point. We have, 
then: 

u x = x X/X u y =yY/Y u z =zZ/Z (3.21) 

where X = dX/dt, Y = dY/dt, Z = dZ/dt. Substituting (3.21) into (3.18), 
we find that the density and pressure profiles in the vapor cloud should be in 
the following form: 
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The constants M and E are proportional to the total mass evaporated and 
total initial energy of the gas cloud, respectively. The normalization constants, 
I\ and J 2 , are given by 


h = Trr [ 7/(7 -i)]r(3/2)/r [ 7/(7 - 1 ) +3/2] 
h = Trr [(27 - 1 )/ (7 -1)] r(3/2)/(7 - i)r [ 7/(7 - 1 ) + 5 / 2 ] 

where T(z) is the gamma-function. 21 

Substitution of (3.21) and (3.22) into the gas-dynamics equations (3.18) 
gives the following equations for the elements of the matrix (3.20): 

d 2 X/dt 2 = -dU/dX d 2 Y/dt 2 = -dU/dY d 2 Z/dt 2 = -dU/dZ (3.23) 

Equations (3.23) are formally identical with the equations of motion of a 
particle in the classical mechanics, where the potential energy is 




x 0 Y 0 z 0 y - 1 

XYZ 


(3.24) 


£(7) = 2 7 J 1 £/( 7 - 1 )I 2 M = (57 - 3 )E/M 

Initial conditions for the set of equations (3.23) can be written as X(0) = Xq, 
dX/dt(t = 0) = 0 and similarly for Y and Z. We assume here that the 
initial kinetic energy of the vapor is much smaller than its thermal energy. 
Introducing dimensionless variables 


£ = X/X$, rj = Y/Yo> C = %/Zo , Vo = Yo/X o, Co = %o/Xo, r = t/to 

where to = Xo(3~ 1 / 2 , and Xo is the largest axis of the initial ellipsoid, we 
transform Equation (3.23) as: 

ft" = m" = CC" = (T?oCo/ftC) 7-1 (3.25) 

£(0) = 1, 77(0) = 770 , <( 0 ) = Co, f (0) = rf(0) = C'(0) = 0 

where £' = d£/dr, etc. Thus, the evolution of the vapor cloud in variables £, 
77, and C is determined by three parameters: 7, rjo, and Co- 

Equations (3.25) have the first integral (in terms of classical mechanics, 
the energy integral). A standard procedure leads to the following relation: 

£' 2 4- i 2 + C' 2 + 2 U = C (3.26) 

where 

t/ft,T?,c) = (7-ir 1 || 7 

and C = 2/(7~ 1). For the special case 7 = 5 / 3 , there exists an additional 
integral of (3.25). This integral has been derived in reference 22 and it can be 
written as 

ft + rj 2 + C 2 = 3r 2 + 1 + 77^ + ft (3.27) 

Using (3.25) and (3.27), it can be shown readily that £(r), r){r), and CM 
are all linear functions of r when r —* 00 . This means that the expansion of 
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the vapor plume becomes inertial as the pressure gradients tend to zero. In 
general, Equations (3.25) have to be solved numerically. Numerical integration 
was performed in references 23 and 24 for a wide range of parameters 770, Co? 
and 7 . The limiting shape of the expanding ellipsoid has been calculated for 
each set of parameters 77 0 , Co? 7- It was shown by Anisimov et al . 22-24 that the 
maximum expansion velocity is reached in the direction of maximum initial 
pressure gradient. It means, in particular, that a disk-shaped vapor cloud 
is transformed into a cigar-shaped one and vice versa. It follows from the 
numerical calculation that the limiting expansion velocity in the z-direction 
is determined mainly by the parameter Co and does not really depend on the 
parameter 770 . In turn, the expansion in the y-direction is governed primarily 
by the parameter 770 . As seen in Figure 3.2, the change of Co by factor of 
30 results only in a few percent change of the ratio Y/X. It appears, then, 
that the plume expansion can be considered as a superposition of two nearly 
independent motions. Both motions become inertial at r > 10 3 - 10 4 . At 
later times, the ratios Y/X and Z/X remain close to their asymptotic values 
at r —► 00 . Some of these asymptotic values are presented in Table 3.4 for 

7 = 7/5. 


CoV 

0.7 

0.5 

0.3 

0.01 

1.219592 

5.486703 

1.472280 

5.821500 

1.965698 

6.461423 

0.03 

1.220410 

3.959977 

1.473599 

4.178862 

1.965206 

4.584814 

0.1 

1.217673 

2.653460 

1.464305 

2.767006 

1.930009 

2.961508 

0.3 

1.205422 

1.732922 

1.430687 

1.773907 

1.834970 

1.834970 


Table 3.1: Values Y/X (upper lines) and Z/X (lower lines) at t —► 00 . 

We can see in Table 3.4 that the expansion of the plume in the z-direction 
is essentially independent of its expansion in the (x, y)-plane. 


3.5 Instability of the vapor-plume expansion 
into a vacuum 

The solution obtained in the previous section is based on the assumption that 
the entropy (per unit mass) is constant over the initial vapor cloud. This 
assumption is physically reasonable. It is clear, however, that the real en¬ 
tropy profile depends on the details of the vaporization process and might be, 
generally, nonuniform. We will show in this section that the stability of the 
vapor plume expansion into a vacuum depends on the initial entropy profile. 
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log T 


Figure 3.2: Evolution of the shape of a vapor cloud during its expansion into 
vacuum. The ratio Y/X is shown as a function of r for different values of 
Vo, Co, and 7 . 

1— no = 0.3, Co = 0.01,7 = 5/3, 2 —no = 0.5, Co = 0.01,7 = 5/3, 

3— 7/0 = 0.3, £0 = 0.3, y — 1.2, 4— 770 =: 0.3, £0 — 0.01, 7 — 1.2, 

5 —770 = 0.7, Co = 0.01, 7 = 5/3 


Instability can arise if the entropy density decreases outward. This fact can be 
used to explain experimentally observed turbulence in vapor plumes. Phys¬ 
ically, the instability mechanism is similar to that responsible for convective 
instability in static stratified media when the entropy decreases in the upward 
direction . 6 

It is easy to see that the group-invariant solution of gas-dynamics equa¬ 
tions (3.22) can be generalized to take into account a nonuniform entropy 
distribution over the initial cloud. One can confirm that solutions similar to 
(3.22) exist when the density and pressure profiles have the form: 

p(x,y,z,t) = p 0 (t)[tp(x,y,z,t)} k 
p(x,y,z,t) = p 0 (t)[xp(x,y,z,t)\ k+1 

During the adiabatic expansion the value: 

5 = p/pi = (po/po)V’ 1_fc(7 “ 1) 

is conserved in each Lagrangian fluid element. The adiabaticity condition 
dS/dt+ u -VS = 0 is satisfied if po = Cp $, where C is a constant. According to 
(3.29), the entropy gradient in the expanding vapor cloud is directed outward 
(inward) if k(y — 1) > 1 (£(7 — 1 ) < 1). For a uniform initial distribution 

of entropy k = 1/(7 — 1 ). 

Notice that the expansion is nonsteady; each fluid element is accelerated 
during the expansion. In a frame connected with a particular fluid element 
this acceleration is equivalent to a gravitational force — 8Md 2 R/dt 2 , 6M 


(3.28) 


(3.29) 
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being the mass of the fluid element. It is known from hydrostatics 6 that the 
mechanical equilibrium of a fluid in a gravitational field is unstable if the 
entropy gradient is parallel to the gravity force. We can expect, therefore, 
that the expanding vapor plume will be unstable when k < 1/(7 - 1), i-e., 
when the entropy of vapor in the external part of the plume is smaller than 
that in the central part of the plume. Destabilization takes place because of 
the buoyancy experienced by fluid elements in the nonuniform inertial force 
field. 

Linear analysis of stability of a spherical cloud expanding into vacuum 
was performed by Book . 25 It was found that the initial perturbations grow 
for all spherical harmonics Z > 0 , provided fc < 1/(7 — 1 ). The relative size of 
the perturbations vary asymptotically like the unperturbed radius. At early 
times, however, when the acceleration is not small, the perturbations grow 
quickly. If Z 1 there is almost exponential growth for t < Ro(M/E J 1 / 2 , the 
number of e-foldings during this time being of the order of Z 1 / 2 . The total 
amplification and the time required to approach the asymptotic state when 
the perturbations “freeze out” both increase with Z. Note that as Z —> oo, 
the instability growth rate diverges. This means that the problem is not 
mathematically well posed. Physically, dissipative phenomena set an upper 
limit on the instability growth rate. 


3.6 Dynamics of vapor expansion 
into an ambient gas 

When the vapor plume produced by laser ablation expands into a gaseous 
medium, the ablated material acts on the surrounding gas as a piston. As a 
result, a shock wave is generated in front of the contact surface. During the 
expansion of the plume the vapor pressure drops and approaches the pressure 
of the ambient atmosphere. The shock wave degenerates into a sound wave, 
and the contact surface stops moving. This stage is known as the stationary 
plume. At this stage the pressure gradient inside the plume is small, and the 
pressure at the contact surface is equal to the ambient pressure, po- In this 
section we will study the effect of various parameters on the shape of the vapor 
plume. We will use the special solution of gas dynamics equations considered 
in Section 3.4. Strictly speaking, this solution is valid only for the expansion 
of gas into vacuum. However, the estimate of the plume shape on the basis of 
this solution yields a reasonable agreement with experiments. 

We consider an adiabatic expansion of ellipsoidal vapor cloud with initial 
dimensions 2 Xq, 2 Iq, and Z 0 . We take the density and pressure profiles in 
the form (3.22). To estimate the shape of the stationary plume, we study the 
motion of the isobaric surface: 


p(x, y, z, t) = Po 


(3.30) 



CHAPTER 3. PROCESSES IN THE VAPOR PLUME 


57 


Using (3.22) we can rewrite (3.30) as 


Po = (£»?<) 7 



y2 ~2\7/(7-l) 

»f C 2 / 


(3.31) 


where p 0 = pol 2 {X 0 Y 0 Z 0 ) / E, x — x/X 0 , y = y/Y 0 , z = z/Z 0 . It can be 
shown that the distance of each point on the surface (3.31) from the origin 
reaches its maximum at some instant t*. We illustrate this in the simplest 
case—a spherical cloud. The radius R = (x 2 4- y 2 4- z 2 ) 1 / 2 of the isobaric 
surface satisfies the equation: 


R 2 (t) = £ 2 (t) [l -po ( 7 ~ 1)/7 £(t) 3(7-1) 


(3.32) 


The equation of motion for £(r) follows from (3.26): 


- 


13(7-1) 


’ (, - **>->) 


1/2 


m = i 


The condition for the existence of a real solution of Equation (3.32) is po < 1. 
Differentiating Equation (3.32), we find 

2 RdR/dT = ft' h, - (3j - l)^ 3 (T- 1 )p 0 (T- 1 )/7 

The derivative dR/dr is positive at r = 0, provided 

Po < [2/(37 - 1)] 7/(7_1) < 1 (3.33) 


If the inequality (3.33) is fulfilled, R(r) reaches its maximum value at r > 0. 
The maximum radius of the isobaric surface is in this case 


Rma X = [3(7 - l)/(3 7 - l)] 1/2 [2/(3 7 - l)] 1 /^-!)^-!^ 

Otherwise, if [2/(37 ~ l)] 7 /^ 7-1 ) < p 0 < 1, the maximum of R{r) is reached 
at r = 0. 

The general three-dimensional case can be investigated in a similar way. 
Equation (3.31) must in this case be solved numerically. The results of cal¬ 
culations and the comparison with experiment can be found in Stangl et al. 26 
Although the model under discussion is extremely simplified, it correctly de¬ 
scribes the dependence of the plume length on the gas pressure. 

The above consideration shows that the contact surface that separates the 
ablated material from the gaseous ambient is decelerated during its motion 
and stops at R = Rmax • In the coordinate frame that moves with the contact 
boundary, this deceleration is equivalent to a gravitational force directed away 
from the ablation products toward the shock-compressed gas. The pressures 
are the same on both sides of the contact surface, while the temperature of the 
shock-compressed ambient gas, under certain conditions, may be higher than 
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that of the adiabatically expanding vapor; therefore, the vapor density may 
be higher than the gas density In this case the contact boundary is unstable. 
This instability is similar to the Rayleigh-Taylor instability (RTI) of static 
equilibrium between two fluids of different densities in a gravitational field. 
The equilibrium is unstable when the density of the upper fluid is higher than 
that of the lower fluid. The RTI is encountered in widely different situations. 
The most well known examples are: supernova explosions, inertial confine¬ 
ment fusion, detonation, cavitation, cumulation, and magnetic field compres¬ 
sion. The instability of the contact boundary between an expanding vapor 
plume and a surrounding gas is quite similar to the instability of the interface 
between the detonation products and the air in a spherical explosion. A linear 
analysis of this instability was performed by Anisimov and Zel’dovich . 27 Near 
the contact boundary, the scale length of the unperturbed flow is h oc c 2 / | g |, 
where c s is the sound speed and g = d 2 R/dt 2 is the acceleration of contact 
surface. Let us examine the evolution of fast-growing perturbations with wave¬ 
length A h. During the perturbation growth time, the acceleration g and 
the density of heavy (pi) and light (p 2 ) fluids all remain practically constant, 
so that the instability growth rate is a slowly varying function of the time. 
Assuming the perturbations of hydrodynamic variables to be proportional to 
the factor P/(cos0) x exp (f 7 dt), where Pi is the Legendre polynomial, and 
considering the limiting case l 1 , we obtain a set of linear equations that 
can be reduced to two equations for the velocity perturbation potential: 


(Pi(r, 6,t) = fi(r) x P/(cos0) x exp (f 7 dt) 
fi + Fi(r)fl + [l 2 /c 2 si 4- l 2 /R 2 } h = 0 


(3.34) 


Here Fi(r ) = dlog(pi)/dr, p*(r) is the unperturbed density profile, i = 1 refers 
to the vapor, and i = 2 refers to the surrounding gas. The quantities p*, c 5 *, 
P, and 7 are slowly varying functions of time. The boundary conditions for 
Equations (3.34) follow from the continuity of the pressure and normal velocity 
component at the (perturbed) contact boundary. Solving Equations (3.34) 
in the semiclassical approximation and imposing the boundary conditions, we 
obtain the dispersion equation for short-wavelength perturbations in the form: 

7 2 = At {l | g/R | +g 2 ( 1 - At 2 ){c~ 2 - c^)/ 4} (3.35) 


where At = (pi - p 2 )/(pi 4- P 2 ) is the Atwood number. The first term in 
the parentheses in Equation (3.35) corresponds to the limiting case of an 
incompressible fluid; the second term takes into account the compressibility 
(to a first approximation in X/h). In the limit l 1, this term does not depend 
on 1. Depending on the ratio of the acoustic velocities, the compressibility 
can have either a stabilizing or destabilizing effect on the contact boundary 
motion. 

According to Equation (3.35), the instability growth rate diverges as l —> 
00 . As in Section 3.5, the instability problem is not mathematically well posed. 
The dissipation and surface tension (which are not taken into account in our 
analysis) lead to the stabilization of short-wavelength perturbations. 
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The RTI of the contact surface is observed experimentally. 14,28 Estimates 
show that under typical experimental conditions, a substantial growth of per¬ 
turbations with l oclO-50 can be expected. These short-wavelength pertur¬ 
bations do not lead to an important asymmetry of the vapor expansion. The 
primary manifestation of a short-wave instability is the mixing of the vapor 
with the ambient gas near the contact surface. The intermixed layer con¬ 
serves the symmetry of the unperturbed flow. The mixing process can be 
described as turbulent diffusion. A detailed analysis of this problem is given 
in reference 29. 
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Chapter 4 


EFFECTS OF ULTRASHORT 
LASER PULSES 


In this chapter we shall concentrate on the interaction of picosecond and sub¬ 
picosecond laser pulses with matter. The problems related to this interaction 
have attracted intense interest since the invention of picosecond lasers. The 
problem of electron relaxation dynamics is of particular interest due to impor¬ 
tant applications of ultrafast laser technology, such as the ultrafast desorption 
of surface adsorbates, generation of ultrashort X-ray pulses, ultrafast melting 
and annealing. 

We will describe the electron and lattice temperature dynamics in a metal 
absorbing an ultrashort laser pulse, the laser-induced electron emission, and 
anomalous emission of short-wavelength radiation from metal surfaces. In 
this analysis, we will consider a metal to be a two-temperature system, with 
different electron and lattice temperatures. We will briefly discuss some of 
the recent experimental studies of hot-electron transport and relaxation in 
metals, which confirm the applicability of this simple model to the study of 
interaction of ultrashort laser pulses with metals. 


4.1 Electron and lattice 

temperature dynamics in a metal 
heated by ultrashort laser pulses 

At energy deposition by an ultrashort visible (or near IR) laser pulse, most 
of the energy is absorbed by “free” electrons through inverse Bremsstrahlung. 
The energy then is thermalized within the electron subsystem. The heating 
of the lattice takes a longer time due to the large ion-to-electron mass ratio. 
As first was indicated in references 1 and 2, the electron gas excited by an 
ultrashort laser pulse can be driven to transient nonequilibrium with the lat¬ 
tice. In this situation, the metal can be considered a two-temperature system 
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and is described by the following set of equations: 1,2 

f CedTe/dt = d/dz(KdT e /dz)-a{T e -Ti) + Q 

1 CidTi/dt = a(Te-Ti) 


Here c e and c* are the heat capacities per unit volume of the electron and 
ion subsystems, respectively; T e and X* are the temperatures of electrons 
and the lattice; and a is the parameter characterizing the energy exchange 
rate between the electron and lattice subsystems. We neglect the phonon 
component of thermal conductivity, since it is much smaller in metals than 
is the electron component. The laser energy deposition Q is taken (as in 
Chapter 1) in the simplest form: 

Q = I(t)( 1 — R)fiexp(—/j,z) (4.2) 

We assume that the laser beam travels along the z-axis and that the surface 
of the metal is located at z = 0. 

At electron temperatures kpT e Ep, where Ep is the Fermi energy, the 
electron thermal conductivity k in the nonequilibrium two-temperature metal 
can be approximated as follows: 3 


K = KoTe/Ti (4.3) 

where ko = constant is the thermal conductivity of metal in equilibrium (T e = 
Ti). The electron heat capacity c e depends on the electron temperature; at 
ksT e Ef , c e « (3T e Ci, (3 = 'K 2 n e k 2 B /2Ep. The energy exchange rate a 
was determined in references 4 and 5. It was shown that a is proportional to 
the electron-phonon coupling constant and can be estimated from the electrical 
conductivity measurements. Estimates for different metals give values of a on 
the order of 10 10 -10 12 W/cm 3 K. 

Equations (4.1) contain two characteristic times, associated with the en¬ 
ergy exchange between electrons and the lattice: r\ = c e /a —the time of 
electron cooling, and 7*2 = c*/a—the time of lattice heating. If the laser pulse 
length ti is shorter than 7 * 2 , the lattice remains cold during the laser pulse. If 
t ^ r 1 , the electron-lattice energy exchange term may be neglected in both 
equations (4.1). In the latter case, Equation (4.1) can be solved easily. Solving 
for the surface temperature, T e (0,£), gives: 

t 1 V 2 

T 0 2 + J I{t — s) exp (ps) erfc(ps) 1 / 2 ds 
0 

where p = p 2 Ko/PTo and To is the initial temperature. At later times, t ti, 
the term c e dT e /dt in the first equation (4.1) can usually be neglected due to 
small electronic heat capacity. At the same time, one also can neglect the 


Te(0,t)= 
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lattice temperature, which is much smaller than the electron temperature, 
Ti <C T e , if t 7*2. The set of equations (4.1) can thus be reduced to: 


d/dz(K,dT e /dz) — aT e + AJ(£)/iexp(—/iz) = 0 

CidTi/dt — aT e = 0 



At low laser intensities, when the parameter K$fi 3 AI / ol 2 T$ < 1, the cooling 
of electrons is mainly due to energy exchange with the lattice. We can ignore 
the heat conduction term in Equation (4.4), and obtain the following simple 
solution for the electron temperature: 

T e (z,t) = I(t)(AjLL/a) exp(-fiz) (4.5) 

In the opposite case, Kofi 3 AI/a 2 To 1, the cooling of electrons and the 
spatial profile of electron temperature are determined by electron thermal 
conduction. The solution of Equation (4.4) may be written in the form: 

T e {z, t) = T e (0, t)(l - z/z 0 ) 2 (4.6) 

where T e (0, t) = [ 3ToA 2 I 2 (t)/2ano J 1 / 3 , and zo = [18«oAZ‘(t)/a 2 To] 1 / 3 . When 
t 7*2, we have T e « Ti, and (4.1) is reduced to the one-temperature heat- 
conduction problem considered in Chapter 1. Approximate analytical solu¬ 
tions of Equations (4.1) in various limiting cases were obtained and analyzed 
in references 6-9. 

Generally, Equations (4.1) can be solved numerically. There are several 
works in which such numerical calculations have been performed to interpret 
experimental results. 10-12 It was found that the model under consideration 
provides a good quantitative description of hot-electron dynamics in metals. 


4.2 Anomalous electron 

and light emission from metals 
heated by ultrashort laser pulses 

Laser-pulse excitation of a metal surface can lead to the emission of electrons 
through two mechanisms: (1) the multiphoton photoelectric effect and (2) 
thermionic emission. Both emission mechanisms have been studied exten¬ 
sively in the nanosecond pulse length range, 13 when the lattice temperature 
Ti was in equilibrium with the electron temperature T e . These experiments 
have demonstrated that the multiphoton photoemission dominates thermionic 
emission for incident fluences less than ~ 0.1 J/cm 2 . Laser fluences above this 
level lead to a significant contribution from thermionic emission. 14,15 In ex¬ 
periments with picosecond laser pulses, multiphoton photoemission has been 
observed 16 primarily at laser fluences below the onset of surface damage. The 
situation changes when excitation times cross into the subpicosecond range. 
The electron subsystem becomes briefly uncoupled from the lattice, enabling 
T e to become larger than T* . 1,2,10 Because the heat capacity of electrons is 
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substantially less than that of the lattice, the degree of nonequilibrium elec¬ 
tron heating can be very high. According to Anisimov et al., 2 the dividing 
line between multiphoton and thermionic emission is estimated to be about 
1 mJ/cm 2 for 100 fs pulses. This estimate is supported by experiments. 10,17 

Multiphoton electron emission is an instantaneous process, i.e., the current 
pulse shows no delay relative to the laser pulse. For thermionic emission, the 
shape of a current pulse is governed by the time dependence of the target 
temperature. In the case of nanosecond pulses, the maximum of the current 
pulse is delayed relative to the laser pulse. Its duration exceeds the laser 
pulse length 77 and the time of delay depends on the thermal properties of 
the target material. With subpicosecond pulses, the electron temperature 
follows almost instantaneously after the laser radiation intensity. The delay 
of thermionic emission pulse relative to the laser pulse is on the order of T\ = 
c e /a, which can be smaller than the laser pulse length. We therefore see that 
with ultrashort laser pulses, a combination of the two fundamental emission 
processes can be observed, whereupon electrons in the high-energy tail of the 
transient thermal distribution further absorb laser quanta and escape from 
the solid. This mechanism, called thermally assisted photoemission, was first 
considered by Anisimov et al. 13 

We now will discuss the effect of transient electron heating on the char¬ 
acteristics of emission current. If the laser intensity is not very high and the 
parameter A = e 2 \E n \ 2 /mhu> 3 1, the perturbation theory can be used 
for the calculation of n-photon photoemission current density. The transition 
probability, as a function of the electron quasi-momentum, has been calcu¬ 
lated by Kantorovich. 18 This probability should be averaged over the Fermi 
distribution of electrons. If the effects associated with heating are ignored, 
the calculation results in the dependence of photoemission current on the laser 
intensity in the form: j n ~ A n ~ J n , where n is the integer part of (ip/hLU + 1), 
and ip is the work function, i.e., n is the minimum number of photons required 
to emit an electron at T e =0. The transient heating produces electrons with 
energies above the Fermi level, so emission can now take place as the result of 
absorption of fewer than n photons. The calculation of total emission current 
in this case was performed by Kantorovich 18 and Ansimov et al. 19 The partial 
currents considered as functions of the electron quasi-momentum and repre¬ 
senting the absorption of n, n — 1, etc. photons were calculated numerically. 
Next, the results were averaged over a Fermi distribution with the tempera¬ 
ture, which depends on laser intensity. As an example of such calculations, in 
Figure 4.1 the emission current from silver irradiated by a neodymium glass 
laser is shown as a function of laser intensity with the electron heating taken 
into account (solid curve). The dashed curve was obtained by ignoring this 
heating ( T e = 0). This line represents the dependence j oc I 5 . Real depen¬ 
dence exhibits a higher degree of nonlinearity, n e ff = dlog j/d\ogI > 5. 

In Figure 4.2 the dependence of the emission current on the laser frequency 
at different intensities is shown. The calculation was performed for a silver 
target. 18 We see that at low intensities, the dependence of the emission current 
on the light frequency is not monotonic. However, at high intensities, the 
electron heating mixes currents of different orders and this smooths out the 




68 


INSTABILITIES IN LASER-MATTER INTERACTION 


-5 



Figure 4.2: Dependence of the emission current on the photon energy for 
different laser intensities (W/cm 2 ): 1—10 9 , 2—4 x 10 9 , 3—7 x 10 10 , 4—10 11 . 
Surface reflectivity R = 0.99. 


The estimates 2,11 and experimental measurements 12,22 show that the elec¬ 
tron temperature in metals irradiated by ultrashort laser pulses can reach 
values about 1 eV, whereas the lattice temperature can remain lower than 
melting point. Under this condition, along with the phenomena of thermoe¬ 
mission and photoemission of electrons discussed above, thermal luminescence 
of the electron gas in the metal also may take place. A significant portion of 
the radiation will be in the visible and UV regions of the spectrum. The delay 
of the luminescence pulse relative to the laser pulse is expected to be on the 
order of T\. 

Experimental observations of the luminescence of metals heated by ultra- 
short pulses was reported in references 8, 23, and 24. Temporal and spectral 
characteristics of the radiation emitted from the silver and tungsten surfaces, 
subjected to 10 and 40 ps laser pulses, were determined. The choice of metals 
used was dictated by the fact that the value of a for W is about two orders 
of magnitude higher than for Ag. Estimates show that in the experiments 
under discussion, electrons in W may be in thermal equilibrium with the lat¬ 
tice, while in Ag a substantial overheating of electrons with respect to the 
lattice can be expected. The measurements made at low laser fluences in the 
560-720 nm spectral range show that the luminescence pulse duration in Ag 
is comparable to that of the laser pulse, while in W the luminescence pulse is 
several times longer. 

The dynamics of the luminescence spectrum were studied at higher laser 
fluences for the silver target in the 510-560 nm spectral range. The most 
intense lines occurring at 521 and 546 nm belong to this range. The fluence 
was 1.0 J/cm 2 , and the laser pulse length was 10 ps. Under these conditions, 
a noticeable evaporation of the target may be observed. The measurements 
show that during the initial 400 ps, the emission spectrum is continuous. 
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Figure 4.3: The luminescence intensities in various spectral bands versus the 
inverse laser fluence (solid lines). The broken lines a and b correspond to 
the l/Q-dependence of the black-body radiation at T = T e and T = T*, 
respectively. The laser pulse duration is 10 ps, the time resolution is 6 ps. 


Later, at t ~0.5-l ns, bands appear, and at t > 2 ns the spectral lines of 
atomic silver are observed. Note that experiments of this type can provide 
important information on the kinetics of vaporization of solids and liquids. 

Figure 4.3 shows the radiation intensity in three different spectral bands 
as a function of the laser fluence. This dependence is presented in the vari¬ 
ables (log I m , 1 /Q ), where I m is the maximum of radiation intensity and Q 
is the laser fluence (J/cm 2 ). The laser pulse duration is 10 ps; the time 
resolution is 6 ps. According to calculations 11 under the conditions of the 
experiment, 24 the maximum electron temperature is a linear function of the 
laser fluence absorbed. This means that the radiation intensity depends ex¬ 
ponentially on 1/T e , which corresponds to the black- (grey-) body radiation. 
The broken lines in Figure 4.3 correspond to black-body radiation intensity 
at temperatures T e (line a) and T* (line b). We see, then, that at low laser 
fluences (Q < Qth ~ 0.6 J/cm 2 ), the luminescence observed corresponds to 
the equilibrium emission of radiation by hot electrons. The emission signal is 
in this case “inertialess”; i.e., it follows the laser pulse shape. However, above 
the threshold fluence Qth, the radiation intensity observed is higher than the 
black-body radiation intensity at temperature T e . Several mechanisms for this 
phenomenon have been suggested in the literature. One of them is a break¬ 
down of the equilibrium Fermi distribution in metals at high intensities. We 
shall discuss this mechanism briefly. 

The distribution function of electrons in metals is described by the Boltz¬ 
mann kinetic equation. The Fermi distribution function is the stationary solu¬ 
tion of this equation, i.e., it causes vanishing of the collision integral. Kats et 
al. 25 showed that for the particles interacting in accordance with the power law 
V(r) oc r~ a (Coulomb forces, Van der Waals forces, etc.), the collision integ- 
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ral vanishes also for power law distributions /(p) oc p s . Similar distributions 
(spectra) are encountered, for example, in the theory of weak turbulence of 
plasma waves. This class of solutions of the kinetic equation is valid in a finite 
region of the momentum space p\ < p < P 2 , and presupposes the existence 
of an energy source and an energy sink outside this region. The properties 
of systems with such distribution functions, especially those properties that 
are sensitive to the presence of high-energy “tail” particles, differ greatly from 
those of equilibrium systems. The analysis carried out by Kats et al. 25 shows 
that in the case of Coulomb interaction, the distribution function of the form 
/(p) oc p 5 / 2 corresponds to constant energy flux in the momentum space and 
is local, i.e., the collision integral converges at both small and large momenta. 

In plasma wave turbulence, the power law spectra are generated when the 
pumping intensity exceeds a certain critical level. Above this level, the spec¬ 
trum is represented by the superposition of the equilibrium (e.g., Rayleigh- 
Jeans) spectrum and the power law spectrum. Similarly, in the case of elec¬ 
trons in a strong radiation field, the power law electron distributions can be 
expected to appear at sufficiently high laser intensities. A detailed analy¬ 
sis of experimental data is necessary to prove (or disprove) this mechanism of 
anomalous emission of radiation from metals heated by ultrashort laser pulses. 
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Chapter 5 


INSTABILITY 
OF SUBLIMATION WAVES 
IN SOLIDS: 

LINEAR THEORY 

5.1 Corrugation instability 

of stationary evaporation waves 
in highly absorbing materials 

In this section we will study the stability of planar sublimation waves in 
strongly absorbing solids. We will show that under certain conditions a pla¬ 
nar, steady-state vaporization front appears to be unstable with respect to 
small perturbations, whose wave vector lies in the front plane. Nonlinear de¬ 
velopment of this instability leads to the formation of nontrivial evaporation 
regimes. 

Although it is more common to investigate evaporation experimentally 
from a liquid phase, we will begin with the study of the solid-to-vapor phase 
transition for several reasons. First, the mechanism leading to the instability 
in the case of solid-to-vapor phase change works also in the case of liquid- 
to-vapor phase change; however, the analysis in the latter case is more com¬ 
plicated. Second, in materials with a rather high saturated vapor pressure 
at the melting point, the layer of liquid near the vaporization front is very 
thin and has only a slight effect on the development of instability. Third, 
a number of materials—e.g., carbon, several plastics and ceramics—do not 
form a liquid phase under conditions typical for laser experiments. A detailed 
linear analysis of instability of the liquid-vapor phase boundary will be given 
in Chapter 7. 

The instability under consideration plays an important role at moderate 
laser intensities, I « 10 5 -10 8 W/cm 2 . As was mentioned in Chapter 1, a sharp 
boundary exists between the vapor and condensed phases at these intensities. 
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Figure 5.1: Schematic showing the vaporization front perturbations. 


The velocity of this boundary, as well as the temperature distribution in the 
condensed phase were calculated in Chapter 1 for the stationary regime of 
evaporation. Now we shall examine the stability of this boundary with respect 
to small perturbations. 1 The situation is shown schematically in Figure 5.1. 
We will introduce small periodic perturbations to the vaporization front posi¬ 
tion and to the temperature field, and will study their evolution over time. As 
in Chapter 1, we will describe the temperature field by the heat conduction 
equation 

cpdT/dt = kAT + Q (5.1) 

where Q = Apl exp[p(Z— z)], and Z = Z(x, y, t) is the position of the solid-gas 
phase boundary. We use below the same notations as in Chapter 1. Boundary 
conditions for Equation (5.1) can be written as follows. As z —> oo , the tem¬ 
perature T approaches To (we will assume here To =0). On the evaporation 
front z = Z(x, y,£), we must impose the continuity condition for the normal 
component of the energy flux. This condition reads: 

kVT = pV n AH (5.2) 

where V n is the normal component of evaporation front velocity, and AH is 
the difference between specific enthalpies of gaseous and condensed phases. 
As in Chapter 1, we can set the enthalpy jump equal to the specific heat of 
vaporization, AH ~ L e //. However, for the nonplanar phase boundary, the 
value of L e f j depends on the local curvature of a phase boundary, and can be 
represented as follows: 

Leff = L v - ( a/p){l/R\ 4- 1 /^ 2 ) (5-3) 

Here, L v is the specific heat of vaporization at the plane boundary, a is the 
surface energy density (surface tension) for the solid-gas boundary, #1 and 
R 2 are the principal radii of curvature of the boundary, and the normal is 
assumed to be directed toward the condensed phase. As we will see later, 
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the dependence of the specific heat of vaporization on the phase boundary 
curvature results in the stabilization of short-wavelength perturbations. 

The evaporation front velocity in Equation (5.2) is determined by the 
kinetics of the transition of atoms from the condensed phase to the gas. When 
a solid evaporates into vacuum, the mass flux is proportional to the saturated 
vapor pressure, and we can use Equation (1.16) for V n . 

We begin our analysis by studying the stability of the stationary evap¬ 
oration wave. The temperature profile in this wave is described by Equa¬ 
tion (1.19). The temperature and the stationary velocity V s of the solid-gas 
boundary can be obtained from the solution of Equations (1.20) and (1.21). 
Introducing new variables: 

£ = MX T) = ny C = m( z - V»t) r = tV?/x 0 = V s /nx 

©(£, V, <» T) = k/xT(x, y, z, t)/AI 8Z(Z, T), t) = n[Z(x, y, t) - V a t] 

we transform Equation (5.1) and boundary condition (5.2) into the following 
form: 


/? 2 d© /dr = pdG/dC + A© + exp(<5Z - C) (5.4) 

V© = pV n L ef f/AI (5.5) 

According to (1.16), the normal component of the evaporation front velocity 
can be written as follows: 

V n = Vo exp {-U eff /e[^6Z(^r),r}} (5.6) 

where U e /f = ML e /f/kB . The planar stationary evaporation wave is de¬ 
scribed by the solution © 5 (C) of Equations (5.4) through (5.6), which depends 
on the variable £ only. Performing simple calculations, we arrive at the result 

©s(C) = Bi exp(-C) 4* B 2 exp(-^C) 8Z = 0 

B x = ((3 - l)" 1 B 2 = -A - Bi/ZJ A = pxuLo/AI 

which coincides with (1.19). As mentioned in Chapter 1, the stationary tem¬ 
perature defined by (1.19) reaches its maximum at finite distance from the 
phase boundary £ = £ m > 0, i.e., the temperature gradient near the vapor¬ 
ization front is positive. Performing simple calculations, we obtain for the 
coordinate of the temperature maximum: 

Cm = 03 - l)" 1 log[l + A (3(0 - 1)] 

The maximum temperature is given by: 

0 sm = /3" 1 [l + A/3(/3-l)] 1 /( 1 - /3) 

It is easy to see that under typical conditions the distance of the maximum 
from the surface is on the order of skin depth, pT l . 
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We will now investigate the stability of the stationary solution (5.7). We 
introduce small perturbations to the temperature and to the phase boundary 
position in the form: 

Q(Z,V,C,t) = © s (C)+<50 <50 = /(C) expiry+ 7T) 

6Z(£,t),t) = E exp(ikr) +-yr) 


Substituting (5.8) into (5.4) through (5.6) and leaving the terms linear in per¬ 
turbations, we arrive at an eigenvalue problem for dimensionless perturbation 
growth rate, 7 . If the real part of 7 is positive, then the solution (5.7) appears 
to be unstable. One should note that the linear growth rate 7 depends only 
upon the modulus of the wave vector, fc. For this reason the dependence of 
6G and 6Z on the variable £ was omitted in (5.7). 

To derive the linear equations for small perturbations, we first should 
project the boundary conditions (5.5) and (5.6) imposed on the perturbed 
evaporation front £ = 6Z(t),t) onto the plane £ = 0 . For this purpose we 
expand the boundary conditions in the power series in terms of the small 
quantity <5Z, and confine ourselves to linear terms. Performing simple calcu¬ 
lations we obtain the following result: 


f" + (3f - (f3 2 y + fc 2 )/ = -Eexp(-C) 


(5.9) 


/'(0) = E( 7 /3 2 + /3AAfc 2 - 0(0, 

/(0) = S(0o/3 7 p _1 — 0q 4- © 0 Afc 2 ) 


as < = 0 (5.10) 


where A = crp/ pLo, ©0 = © 5 ( 0 ), ©0 = ©^(0), etc. For a rough estimate 
of the order of magnitude of A, note that if we substitute in Equation (5.3) 
R~k~ x ~d,d being the interatomic distance, then both terms on the right 
side on this equation must become on the same order. We immediately obtain 
the estimate A ~ fid 1. 

Now, we can solve Equation (5.9) with the boundary conditions (5.10) 
to obtain a dispersion equation which defines the dependence 7 (fc ). 1 This 
equation can be written conveniently in parametric form as: 


_ ct — f3 ©q — (a — (3 4-1) 1 — aA(a©o — © 0 ) 

7 /3 a0 o (0oP -1 - A P) + (A/3(l - A/3) ( 5 . 11 ) 

fc 2 = a(a —/3) —/3 2 7 Rea > 0 


where a is a parameter and p = XcM/ks . 

We begin our analysis of the dispersion equation (5.11) with the case of 
real values of a, which correspond to real values of 7 . Consider, first, the 
denominator of Equation (5.11). Estimate the magnitude of the product A/3. 
Using the approximation A ~ fid and changing to dimensional variables, we 
find that A/3 ~ V s d/x • But x/K is the spatial scale of the temperature field 
in the stationary evaporation wave (see (5.7)). For a macroscopic description 
of the problem to make sense, this quantity should be much greater than the 
interatomic distance. We see, hence, that A/3 1, and 1 — A/3 > 0 . One can 
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show, also, that the product A/3 <C ©o p 1 • Indeed, the inequality ©o pX0 
can be rewritten as 


©o P 1 exp(p© 0 J ) » XVq/hx ~ Vo d/x (5.12) 

A typical value of d is ~ 10“ 8 cm, Cq ~ 10 5 cm/s, and x ls not less than 10 -2 
cm 2 /s for condensed material. We see that the right side of Equation (5.12) 
is small compared with unity, which makes the inequality obvious. Thus, we 
see that the denominator in Equation (5.11) is positive. 

We will now investigate the numerator of Equation (5.11). First, we see 
that, when a = (3, 7 = k = 0. Further, since ©q = X/3 < 1, we have at 
0 < (a - (3) < 1 

©0 - (a - 0 + l)" 1 « \(3 - 1 < 0 (5.13) 

It follows that 7(fc) < 0 for small fc, i.e., the long wavelength perturbations 
decay with time. As a increases, the magnitude of the term (a — /3 - fl)" 1 
decreases and when 

a — a\ = /? — 1 4* l/©o (5-14) 

the left side of Equation (5.13) vanishes. If the latter term in the numerator 
aA(a© 0 — ©0) is small at a = a 1? in comparison with ©q, then the growth 
rate 7 changes sign at a ~ aq. Positive value of 7(fc) means that a plane 
stationary evaporation wave becomes unstable. 

At large a, the term aA(a©o — ©0) is large in absolute value. Therefore, 
the growth rate 7 again becomes negative, if a > c*2. A reasonable estimate 
for the magnitude of a .2 is: 


a 2 ^ (©o/A© 0 ) 1/2 (5.15) 

If A is not too small and aiA(ai©o — ©0) > ©o> stabilizing effect of the 
last term in the numerator of (5.11) asserts itself at a ~ aq, and 7(fc) may 
remain negative for all finite values of k. 

Using the estimate for the magnitude of A, we may readily verify that Acq 
is always small compared with unity. Hence, we may conclude that the plane 
stationary evaporation wave is stable when a 2 A©o > ©0? and unstable when 
a 2 A©o -C ©o- In the latter case, small perturbations of plane vaporization 
front grow with time if their wave numbers lie in some finite interval k\ < k < 
fc 2 , where k\ ~ ai(ai - /?), fc 2 — <*2(^*2 - (3 ). 

Returning to the solution of the stationary problem, one can show easily 
that ©0 decreases monotonically as the laser intensity increases. This means 
that the stabilizing term a 2 A©o related to the surface energy also decreases 
with laser intensity increasing. Therefore, a threshold in laser intensity must 
exist, above which a plane stationary evaporation wave becomes unstable. The 
stability of this wave at laser intensities below the threshold is due entirely 
to the contribution of the surface energy to the latent heat of vaporization. 
If this contribution is neglected, short wavelength perturbations, which make 
the plane stationary evaporation wave unstable, always exist. 
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Figure 5.2: The dependence of the instability growth rate 7 on the perturba¬ 
tion wave number k at different laser radiation intensities: 

1— I = 0.35/*; 2—1 = /*; 3 —I = 3.33/*; /* = 3.33 x 10 ~ 3 pV 0 L v , A = 10" 3 . 


At laser intensities well above the instability threshold, upper and lower 
boundaries of the instability interval in k may be determined by Equations 
(5.14) and (5.15). Within this interval the dispersion equation (5.11) may be 
simplified as: 

7 (k) = (p//30o)[0' o - (k + l) -1 - fc 2 A0 o ] 

From the last equation, one can see that a maximum of the instability growth 
rate is reached at k = k m , where: 


and is equal to 


k m ~ (2A0 O )- 1/3 


7m = (p//30q) 



1.9(A0 O ) 1/3 


(5.16) 

(5.17) 


Detailed numerical investigations of the dispersion equation (5.11), made with¬ 
out simplifying assumptions, confirms the simple analysis described above. In 
Figure 5.2 the dispersion curves 7(fc) are shown corresponding to three differ¬ 
ent values of the laser intensity: subcritical / < /* (curve 1), critical 1 = 1* 
(curve 2), and supercritical / > /* (curve 3). In the last case, instability 
growth rate is positive in a finite interval of wave numbers ki < k < fc 2 . 

We now turn to Equation (5.17) to find the maximum growth rate 7 m . 
Because the negative term in (5.17), proportional to A 1 / 3 , is significantly 
less than the positive term in the region of instability, we have the following 
estimate: 

7m«P©o// 00 O ( 5 - 18 ) 

Transforming (5.18) to dimensional variables and using Equation (5.6), we 
obtain: 


7m * (dV s /dT s )T' s = dV s /dz 


(5.19) 
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We see from (5.19) that the instability growth rate is positive if a forward 
displacement of the evaporation front causes its velocity to increase. 

We have studied the branch of dispersion equation (5.11), which corre¬ 
sponds to real 7 and real fc. This branch describes a nonoscillatory (aperiodic) 
growth of small initial perturbations. One can show, using Equation (5.11), 
that complex branches of the dispersion equation lie in a nonphysical domain 
for the case of strongly absorbing materials. For complex 7 , the condition 
that the temperature perturbations must decrease at infinity, Rea > 0, is not 
satisfied. This means that only the aperiodic growth of small perturbations is 
possible in the case of sublimation of strongly absorbing solids. In Chapter 4 
and in Section 5.3 of this chapter, we shall show that an oscillatory growth 
of small perturbations is possible for laser evaporation of liquids and laser 
sublimation of dielectrics. 


5.2 Corrugation instability 

of nonstationary evaporation waves 
in highly absorbing materials 

As it was shown in Section 1.3 (see also references 3-4), for highly absorbing 
materials irradiated by a laser of constant intensity, the stationary regime of 
evaporation is reached within the time on the order of To ~ x/K 2 - Transform¬ 
ing (5.18) to dimensional variables, we obtain: 

7m » ( V?/x)(M/ckg ) [Lo/T s (0)} 2 » To 1 

Thus, when the laser intensity is significantly above the threshold and Equa¬ 
tion (5.18) holds, the growing perturbations could destroy the plane phase 
boundary before the stationary mode of evaporation is established. In this re¬ 
gard, a question arises about how perturbations of a plane vaporization front 
approaching a stationary mode evolve. To examine this question, we shall 
now investigate the development of small perturbations of plane nonstationary 
evaporation front and corresponding temperature field. We shall assume that 
the laser radiation is switched on at time t = 0 , and has a constant intensity 
J, which considerably exceeds the instability threshold of a plane stationary 
wave. From the previous discussion we may expect that there are two different 
types of evaporation front motion: “slow” one-dimensional motion, associated 
with a change of characteristics of the plane vaporization wave, and “fast” mo¬ 
tion, associated with the development of non-one-dimensional perturbations. 
This enables us to investigate the problem in an adiabatic approximation, i.e., 
to study the stability of one-dimensional motion with respect to perturbations 
of the form: 


m m 

T 

ikrj 4 - / 7 (r ; ) dr f 

6Z = Eexp 

r 

ikr) -b / 7 (r ; ) dr f 

J 

0 


J 

0 


<50 = /(C) exp 
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assuming that 7(7*) is a slowly varying function of time: 

\d^/dr\ < 7 2 

The last condition means that the change in the instability growth rate in the 
characteristic time of the instability’s development should be small. 

It is convenient to consider the instability problem in a coordinate frame 
connected with an unperturbed evaporation front. Since in the present case 
the front velocity changes with time, we must slightly modify our definition 
of dimensionless variables. We assume now 


C = m 

t 

z- f V(t') dt' 

6Z = (i 

t 

Z(y, t) - [ V(t') dt' 


J 

0 


J 

0 


The quantity x/V? is no longer a characteristic constant, so that the dimen¬ 
sionless time r is now chosen in the form 

r = 2 

The remaining notations agree with those used previously in this chapter. 

The analysis, which is similar to that described above, leads to the disper¬ 
sion equation in the form: 

©0 + (a - P) [A/3 - (a - (3 + l )" 1 - aA(a © 0 + A f3)] 

7 = 13 a0 o (0 O p - 1 - A/5) + A/3(l - A/3) (5.20) 

k 2 = a(a — (3) — 7 Rea > 0 

Here ©0 = ©o(t) is the temperature on the surface of the condensed phase, 
the dot indicates a derivative with respect to time, and (3 = /3(t) = fiV(t)/Xi 
where V(t) is the instantaneous velocity of the unperturbed plane evaporation 
front. The asymptotic values of the front temperature and velocity as r —► 00 
(corresponding to the stationary evaporation wave) will be marked by the 
additional subscript s, i.e., they will be denoted as ©o 5 and f3 s , respectively. 
Note that 7 in Equations (5.11) and (5.20) differs by the factor (3 2 due to the 
difference in the definition of dimensionless time. 

Despite the formal resemblance between the dispersion equations (5.11) 
and ( 5 . 20 ), the spectrum of perturbations, which develop on the background 
of a nonstationary wave differs substantially from the corresponding spectrum 
in the stationary problem. The difference is because (3(t) in (5.20) varies with 
time: (3(0) = 0 and /3(t) —► /3 S as r —► 00 . Since A is independent of r, 
the second term in the numerator of Equation (5.20) is negative at small 
r. This means that 7 (k) decreases monotonically as k increases. Thus, the 
plane nonstationary evaporation wave should be stable at significantly small 
r for any value of laser radiation intensity. Since f3(r) increases with time, the 
contribution of the term A (3 in the numerator of Equation (5.20) increases and 
at a certain time r = To the growth rate 7 (k) becomes positive, first at some 
point k = ko ^ 0 . Afterward, the instability region expands to a finite interval 
and the value of 7 continues to grow. As an example, a family of dispersion 
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Figure 5.3: Dispersion curves for a nonstationary evaporation wave. Laser 
intensity / = 10~ 4 pVoL Vl A = 10“ 5 . 

1 — r = lO” 7 , (3 = 3.4 x 1(T 5 , T = -0.11, 

2— r = r 0 = 1.7 x 10~ 7 , (3 = 4.6 x 10" 5 , T = -0.15, 

3— r = 2.5r 0 , (3 = 6.2 x 10" 5 , F = 0.35, 

4— r = 3.8r 0 , p = 6.6x 10" 5 , F = 1.10, 

5— stationary regime, f3 = (3 S = 7.5 x 10“ 5 . 


curves is shown in Figure 5.3 for different times (ri < tq, T2 = r 0 , T3 > r 0 ) 
along with the curve corresponding to the stationary evaporation wave r —> 00 . 
Note that the function 7 (&) reaches its maximum when a (3. Hence, for 
the position of the maximum fc m , Equation (5.16) holds km « (2A©o(r))“ 1 / 3 . 
This means that the maximum of the growth rate is displaced with time. The 
growth of perturbations with time may be characterized by the integral 



For curves 1 , 2 , and 3 in Figure 5.3, the values of F are - 0 . 11 , -0.15, and 0.35, 
respectively. We may estimate the value of To from the condition: 

max[7(fc,T 0 )] = 0 
k 

Using Equation (5.17), we obtain the following equation for tq: 


A0(t o ) « [A0 o(t o )] 1/3 


( 5 . 21 ) 
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To determine r 0 explicitly, we should know the dependence P(r). This depen¬ 
dence was found using a numerical solution of the heat conduction problem 
(1.17) and (1.18) (see references 5-6). To obtain an order-of-magnitude esti¬ 
mate of r 0 , we use the simplest analytical approximation for /?(t), which is in 
reasonable agreement with the numerical calculations 

P(t) = P s exp(-l/r/?J) (5.22) 

Substituting Equation (5.22) into (5.21) and assuming that X/3 « 1 , we obtain 
for To: 

3 

T ° ~ PI log(A0 Os ) 

The characteristic time t* of the destruction of the plane nonstationary evap¬ 
oration wave can be determined from the condition 


T(t*) = max 
1 k 



1 


It is easy to see that due to the rapid increase in 7 (t, k) when t > t 0 , t* must 
be close to t 0 . Thus, despite the fact that the laser intensity I considerably 
exceeds the threshold intensity /*, the time required for destruction of the 
plane nonstationary wave is of the same order of magnitude as the time to 
form the stationary evaporation wave. 


5.3 Instability of plane stationary 
evaporation waves in dielectrics 

In this section we shall study the instability of laser vaporization of mate¬ 
rials whose absorption coefficient is strongly temperature-dependent . 7 Laser 
breakdown and vaporization of nonlinearly absorbing materials have been con¬ 
sidered in Chapter 2 . It was established that to maintain the vaporization pro¬ 
cess, a laser intensity well below the optical breakdown threshold is required. 
To initiate a vaporization wave at this intensity, an absorbing “seed” is neces¬ 
sary, which ensures a local temperature increase and creation of an absorbing 
layer with a sufficiently high free electron concentration. A qualitatively simi¬ 
lar situation occurs in gases , 8 in which the ionization waves supported by laser 
radiation may propagate at laser intensities considerably below the breakdown 
threshold of a cold gas. 

Now, we shall recall some of the properties of vaporization waves in trans¬ 
parent dielectrics. First, the calculations 9,10 show that the optical absorption 
in such waves is small, A = I a bs/Io ^ 1? decreasing with increasing laser 
intensity. In other words, a considerable fraction of the incident light passes 
through the wave not being absorbed. Second, the geometric thickness of the 
absorbing layer in the dielectric is much less than the characteristic scale of 
the temperature distribution. To simplify the instability analysis, we shall 
neglect a small change in the light intensity inside the absorbing layer. The 
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calculation of the evaporation wave parameters is reduced, in this case, to 
determining the temperature field inside the dielectric. 

As mentioned in Chapter 2, the evaporation waves were observed experi¬ 
mentally on both the incident and exit surfaces of the target. In the limiting 
case of very small absorption, the both experimental situations are described 
by the same equations. 

Consider a uniform laser beam with constant intensity, J 0 , travelling 
through a dielectric along the z-axis. The temperature of the dielectric satis¬ 
fies the heat conduction equation: 


cpdT/dt = div (kS7T) + /i(T)J 0 (5.23) 

with the boundary conditions 


kVT = pL v V n (T) and V n = V 0 exp (-U/T) (5.24) 

on the vaporization front z = Z(y,t), and 


T —> 0 at z —> oo 


(5.25) 


As in the previous section, we assume that the latent heat of vaporization 
depends on the curvature of the boundary and is given by L v — Lo + ° / P-R? 
where R is the local radius of curvature of the boundary. We consider vapor¬ 
ization into a vacuum, and the vapor is assumed to be transparent for the 
laser light. For the temperature dependence /i(T) of the absorption coefficient 
of the dielectric, the same expression is used as in Chapter 2. This choice is 
justified because the experimental dependence of the stationary vaporization 
wave parameters on the laser light intensity is in satisfactory agreement with 
calculations described in Chapter 2 (for a constant thermal conductivity, see 
references 9 and 10). 

The boundary value problem in Equations (5.23) through (5.25) has a par¬ 
ticular solution describing a planar stationary evaporation wave: T(y,z,t) = 
r 5 «), where C = z — V s t. To consider the stability of this solution, let us 
introduce small perturbations of the phase boundary and temperature field in 
the form: 


6Z(y,t) = Z(y,t)-V s t = 3 exp(yt + iky) 

6T(£,y,t) = T(z,y,t) -T s (£) = f(Qexp('yt +iky) 


Substituting (5.26) into (5.23) and confining ourselves to linear terms in the 
small quantity /(C)? we find that 


Xf" + v.f + f[xU(0 - 7 - X* 2 ] = 0 (5.27) 

Here 

U(0 = (I 0 /KW/dT\ T=Ta 


and x = K / C P 1 s the thermal diffusivity. As in the previous section, we have 
to project the boundary condition (5.24) related to the vaporization front 
z = Z(y,t) onto the C = 0 plane. Having expanded Equation (5.24) in the 
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power series in terms of the small perturbations 6Z and <5T, we obtain in the 
linear approximation: 


where 


/'(0) + a/(0) = 0 

V. V? + IX + A(fcx) 2 
x V 2 - me[ej X + A(fc X ) 2 ] 


(5.28) 

(5.29) 


A - aV s /pL 0 x e = kBT s (0)/L o m = cM/ks 


ks and M are the Boltzmann constant and molecular mass, respectively. We 
note that unlike the case of linear absorption considered in Section 5.1 the 
perturbation of phase boundary 6Z does not appear in the equation for the 
temperature perturbation, thus the analysis of stability is somewhat simpli¬ 
fied. 

The instability growth rate 7 (k) is the eigenvalue of the problem (5.27) 
and (5.28). It is convenient to transform Equation (5.27) into the Schrodinger 
equation. Setting f(Q = </>(£) exp(— (V s /2x) in (5.27), we find 


4>" + (e-U)<t> = 0 


(5.30) 


where e = k 2 4 7 /x — {V s /2x) 2 - Boundary condition (5.28) then becomes 


</>'(0) -f w<t>( 0) = 0 a ) — ol- V s /2x 


The temperature derivative of the absorption coefficient is analogous to the 
potential term in the Schrodinger equation (5.30). The basic temperature 
dependence of /i(T) is described by the factor exp(—E/T), where E T. 
It is easily seen that because of the strong temperature dependence of /i(T), 
the function U(Q in (5.30) is nonzero only in a narrow interval close to the 
maximum of the unperturbed temperature T S (Q. The determination of the 
eigenvalue 7 (k) is reduced, in this case, to solving a standard problem (well 
known in quantum mechanics) on the energy level in a short-range potential 
field (a slight difference is in the different boundary condition). Straightfor¬ 
ward calculations lead to the following system of equations for the instability 
growth rate 7 (fc): 


2p(p - u) = J\p - u) 4- {p 4- w) exp(-2pCo)] 

P 2 = (V S /2 X ) 2 4 7 /x 4 k 2 u> = a- V s /2x Rep > V s /2x 


where 00 

j = (h/x) j ATs)d<; 

0 

Co is the coordinate of the maximum of the function T S (Q, and a is given by 
(5.29). The set of equations (5.31) can be reduced to a single equation for the 
parameter p, 

(2p - J)SL 4 JS+ exp(-2pCo) = 0 (5.32) 

where 

S± = (p± V s /2xf ~ (1 - A )k 2 + me 2 ( 1/2 T PX/K)-? 0 
So = P 2 ~k 2 { 1 - A/e) - 04/2*) 2 
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kx/V 8 


Figure 5.4: The regions of oscillatory (shaded) and aperiodic perturbation 
growth. 1—Co V a /x = 1, 2—Co Vi/* = 0.1. 


In terms of the solution of Equation (5.32), the instability growth rate is 
given by 

7 = X(P 2 ~ k 2 ) - V, 2 /4x- 

The inequalities A <£ me 2 1 and Co being on the order of, or smaller than, 
x/I 4 , which follow from an analysis of the unperturbed solution 9,10 are used 
to solve Equation (5.32). 

It is easy to see that the range of parameters, within which the growth rate 
7 (k) has a nonzero imaginary part (oscillating solutions), coincides with the 
range of parameters, in which Equation (5.32) has complex roots. This region 
is shown in the (J, fc)-plane in Figure 5.4. The location of the boundaries of 
this region depends on the parameter Co¬ 
in the study of instability, we are interested in modes with Re 7 > 0, which 
correspond to the solutions of Equation (5.32) with a sufficiently large real 
part, Rep > Imp, &, and V s /x- Consider, first, the limiting case Re(pCo) 

1 , when the second term in (5.32) gives a small correction. Neglecting this 
term, we obtain two equations: 

2p - J = 0 = 0 

One of the branches of the solution is, thus, given by 

p 1 = (J/2)(l + il 1 ) (5.33) 

where 

Ri « (-1 + me 2 xJ/V a ) exp(-J^o) 

The solution (5.33) exists when Re(pCo) 1 and R\ 1 . The last condition 
is not fulfilled if pi is close to the solution of the equation S _ = 0. The growth 
rate at p = pi is given by 

7(*) = x(«J 2 /4 - k 2 ) - V 2 /4 x 


( 5 . 34 ) 
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Perturbations with k = 0 have the highest growth rate; as k increases, the 
perturbations become stabilized due to heat conduction. The mode (5.34) is 
similar to a thermal explosion or optical breakdown initiated by a microinclu¬ 
sion (see Chapter 2). 

We note that the wave number k does not appear in the approximate 
equation (5.33). However, it is easily seen directly from Equation (5.32) that 
the solution (5.33) is inapplicable in the exponentially narrow region defined 
by the inequality 

|07 - VJ X ? ~ 4(1 - A )k 2 + me^xJ/Vs + 1) 
x [J 2 - 4(1 - A/e)k 2 - (V s /x) 2 ]| < (4JV s /x)exp(-J( 0 ) 

In this region the growth rate 7 (k) assumes complex values. 

The other branch of the dispersion equation for 0 1 describes unstable 

modes similar to the ones studied in Section 5.1 (corrugation modes). One can 
find from Equation (5.32), P 2 = P 2 o(l 4- # 2 )? where P 20 satisfies the equation 
5 _(p 2 o) = 0- For solving this equation, it is convenient to transform 5_ into 
the following form: 

5_ = (1 -f A) 


V. 1 

2 

- k 2 

A(14A/c) 


V S A 1 

[ P 2x(l + A)J 

— rv 

(l + A) J 


[2 X (1 + A)J 


where A = me 2 ( 1/2 4 - px/V s )- Taking into account the inequality me 2 1 
the solution of the equation S- = 0 can be written as: 

P 20 = Vs/ 2x(l 4 A) 4 k[ 1 — A(1 4 A/c)/(l 4 A)] 1 / 2 4- # 20 ? (5.35) 

where the correction 

R20 ~ [F r 5 A/2x(l 4 A)] 2 /[p — P20 — V s / 2 x{l 4 A)] 

is always small. In limiting cases of small and large wave numbers k one can 
obtain from (5.35) 


P 20 


V s /2x 4 k(l — A ) 1//2 
k( 1 - A/e) 1 / 2 


for k V s /me 2 x 
for k V s /me 2 x 


The condition Rep^o 1 is fulfilled for this branch if k £0 1. The correction 

R 2 is approximately: 

4JV s [V a 2 -me(l-e)AxW} 

2 ~ (1 - 2 p 20 )xk(V s + 2 me 2 X k) 2 P( 

Neglecting the corrections R 20 and # 2 ? as well as A < 1 , we obtain from 
(5.35) the following expression for 7 (k) 




kV s 

(l + A) 


+ A(2 + A) 


2 x(l + A)J 



(1 + A/e) 
(l + A) 


(5.36) 
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with A « me 2 xk/V s . At small and large values of k we have 



kV s - A X k 2 

/me 2 x - Axk 2 /e 


for 1/Co < k < V 0 /me 2 x, 
for k V s /me 2 x- 


The dispersion equation (5.36) is similar to that obtained in the case of 
strongly absorbing materials (see Section 5.1). For this branch, 7 (k) has a 
maximum at finite k ^ 0 . As in sections 5.1 and 5.2, the dependence of 
the heat of vaporization on the curvature of the phase boundary results in 
stabilization of the short-wavelength perturbations. 

Approximate formulas (5.33), (5.34), and (5.36) are not valid in the small 
range of k « J/2, where the two branches of dispersion equation are close to 
one another. Consider this region, first, in the case RepCo 1. The correction 
proportional to exp(— 2pCo) is important in this case only in the narrow range 
of wave numbers \k — &o| V s /x, where ko is determined from the equation 
P 2 o(ko) = J/ 2. In this range, an approximate solution of Equation (5.32) can 
be obtained in the form: 


P ^ (2^20 - J/ 4) 



2 P 20 - Jy JV S [AV 2 - 7716(1 - c)A x 2 J 2 } 


%X[V 2 4- me 2 x J ] 2 


exp(-JCo) 


It is easily seen that p and 7 are complex when k is close to ko. When 
J Co 1? we have: Re 7 Im 7 , i.e., the steady-state solution is unstable 
and the instability is oscillatory. The dispersion dependence 7 (k) is shown 
in Figure 5.5. For both branches of the dispersion equation (5.32) considered 
above, Re 7 is positive in certain ranges of k. The maximum growth rate has 
the thermal explosion mode at k = 0 . 



Figure 5.5: The dispersion dependence 7 (k). m = 3, e = 0 . 1 , A = 10 5 
J = 18 Vs/x- 1— Co = x/V a , 2— Co = 0.1x/Vs. 
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We now consider the case when J( o 1. Solutions of the thermal- 

explosion type, described by Equation (5.34), no longer exist. Aperiodic 
modes of the type (5.36) now are present only for k( o 1. No solutions 

with Re(pCo) 1 exist when o 1 . 

Consider the solutions satisfying the condition |p(o| 1. Using the series 

expansion for the exponential function, we find (for o 1 , J Co ^ 1 ): 

P = (1 - 2JCo)K/2* ± [A : 2 - V,J{ 1 + J<6)/x] 1/2 (5-37) 

For small k , Equation (5.37) gives complex p, which correspond to oscillating 
perturbations with growth rate 

7 (*) = ~V.J ± V s (k 2 - V a J/ X ) 1/2 (5.38) 

Solutions described by (5.37) with \pCo\ ^ 1 also can be constructed for the 
values *7Co 1, &Co 1- Note that perturbations with &Co > 1 and \pCo\ < 1 
clearly have Re 7 < 0 , i.e., they decay. Perturbations, for which |Imp| 
Rep, also are damped and not considered here. 


5.4 Stability analysis of plane stationary 
ablation waves in polymers 

In Section 2.3, we studied the model of UV-laser ablation of organic polymers. 
The model describes the possibility of “cold” ablation, i.e., it yields realistic 
ablation rates at moderate surface temperatures, which are in agreement with 
many experimental observations. It was shown in Section 2.3 that the tem¬ 
perature distribution in the stationary ablation wave has a maximum behind 
the ablation front. Similar maxima appear also in laser evaporation of met¬ 
als and transparent dielectrics and can lead to the corrugation instability of 
vaporization front. In this section the stability of the plane stationary abla¬ 
tion wave will be studied. It will be demonstrated that in polymer ablation, 
excited species stabilize the ablation front when the thermal relaxation times 
are of the order of vp > 10“ 10 s. 

We shall describe the ablation process by the equations for the number 
density of excited particles iV*, the laser intensity, and the temperature, T. 
These equations can be written as: 

dN*/dt = VdN*/dz + ( aI/hv)(N 0 - 2 N*) - N*/t t (5.39) 

dl/dz = -aI(N 0 -2N*) (5.40) 

dT/dt = VdT/dz + x(d 2 T/dz 2 + d?T/dy 2 ) + (x/*)Q (5-41) 

Here we use the same notations and scaling values as in Chapter 2. 

We now investigate the stability of the stationary solutions iV*(z), / 5 (z), 
T s (z) of Equations (5.39) through (5.41). We consider the perturbations of 
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the ablation front 6Z(y,t) shown in Figure 5.1. The boundary conditions on 
the perturbed ablation front are written in the form: 


«V„T \ 2mSZ(y<t) = p [zA/TlT + (1 - x)AHU}\ 2 _ szm (5.42) 

where ftV n T is the normal flux component for the perturbed ablation front, 

x(z,y,t)=N*(z,y,t)/N 0 , U* = V 0 * exp[-E*/T6Z(y,t)}, 

II = Vo exp[— E/T6Z(y, £)], and A if* and A H are the enthalpies of vaporiza¬ 
tion of ground-state and excited-state species, respectively. The velocity of 
the perturbed ablation front is given by: 


V = V s + dsz/dt = [xir + (1 - x)n]\ z=SZ(yt) 

Note that the activation energies E and E * are assumed now to be dependent 
on the local curvature of ablation front. We shall use the same approximation 
as in Sections 5.1 through 5.3: E = Eo — oM/pksR , and similar expression for 
E*, where a is the surface tension (surface energy density), M is the molecular 
mass, and R is the local radius of curvature of the perturbed ablation front. 
We assume for simplicity that a does not depend on the concentration of 
excited species. 

The perturbed quantities can be written as 


6Z(y, t) = Z\ exp (7 1 + iky) 

N*(z,y,t) = N*(z) + Ni(z) exp^t -f iky) 

T(z,y,t) = T s (z) + Ti(z) exp(yt + iky) 
I(z,y,t) = I s (z) + Ji(z)exp( 7 * + iky) 


(5.43) 


Substituting (5.43) into (5.39) through (5.41) and performing linearization 
with respect to the small perturbations Zi, N\, Ti, and ii, we obtain a set of 
ordinary differential equations. The boundary conditions for these equations 
can be obtained from the exact boundary condition on the ablation front 
z = 6Z(y,t) employing the standard procedure of expansion of (5.42) in series 
in the small quantity £Z. The instability growth rate 7 (k) is an eigenvalue 
of the resulting boundary value problem. The lengthy but straightforward 
calculations of the dispersion equation 7 (k) were performed by Luk’yanchuk 
et al . 11 We shall now consider the results of these calculations. 

The dependence 7 (k) is shown in Figure 5.6 for J 0 = I b = 10 7 W/cm 2 
and for various values of r^. The time unit used for the normalization of 7 
and tt is to = 10 ~ 8 s. The parameters employed in the calculation of 7 (k) are 
typical for organic polymers. The figure shows that the dependence 7 (k) for 
the polymer ablation is qualitatively similar to that obtained in Section 5.1 for 
the evaporation of metals. In both cases, a finite region of wave numbers exists 
where 7 (k) is positive and the perturbations increase with time. However, the 
instability appears only when the relaxation process is rather fast. For the 
parameters employed in Figure 5.6 the critical value of relaxation time is tt = 
20 ps. When the relaxation time is longer, the ablation front becomes stable. 
The critical relaxation time depends on the laser intensity. In Figure 5.7 the 
region of parameters Io/h and rx/to is shown where the instability can occur. 
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Figure 5.6: The dispersion dependence for different values of relaxation time 
(in units 10 -8 s). The laser intensity Iq = h = 10 7 W/cm 2 . 



Io/Ib 


Figure 5.7: The region of instability (shaded) in the (Jo,7Y)-plane. 
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It is seen that for relaxation times longer than the critical time Tj.^4x10 -11 s 
the ablation front is stable at any laser intensity. For different polymers the 
critical relaxation time is between 10“ 10 and 10“ 11 s. So we consider tt > 
10“ 10 s as a sufficient condition for stability of the ablation front. 

Now we shall discuss the physical mechanism of stabilization. Returning to 
Equation (5.19), we see that maximum instability growth rate is 7 m ~ dV s /dz. 
It is positive when the temperature grows with z. Considering the ablation 
model proposed in reference 11 , the ablation velocity is a function of both 
temperature and concentration of the excited species. Thus, the instability 
criterion is: 

dV s /dz\ z ^ = (dV s / dT)(dT s / dz) + {dV s /dN*){dN;/dz) > 0 (5.44) 

If E* <C E the derivative dV s /dN* > 0, and the second term in Equation 
(5.44) becomes negative because dN*/dz < 0 (see Chapter 2, Figure 2.10). 
Thus, this term that originates from the activated desorption of excited species, 
stabilizes the ablation front. 

From a physical point of view, the suppression of the ablation front insta¬ 
bilities can be understood as follows: 11 with increasing relaxation times the 
contribution of excited species to the total ablation rate increases. Hence, 
the effective activation energy and, thereby, the value of dV s /dT decreases. 
Additionally, as is seen from (5.44), the dependence of the ablation velocity 
on the concentration of excited species leads to stabilization of the ablation 
front. 
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Chapter 6 


NONLINEAR EVOLUTION 
OF INSTABILITY 
IN LASER EVAPORATION 


6.1 Slightly supercritical structures 
in laser sublimation waves 

It was shown in Chapter 5 that, under certain conditions, a plane stationary 
evaporation front can be unstable. The instability occurs when the laser radia¬ 
tion intensity exceeds a threshold which depends on thermophysical properties 
of the material evaporated. At supercritical laser intensities, a finite interval 
of wavenumbers &i 7 < k < & 2 7 exists in which the perturbations of plane front 
increase exponentially with time. In the weak supercritical regime, i.e., when 
the laser intensity is above, but close to the threshold, the band of unstable 
modes is narrow and located near the critical wavenumber fc m , as it is shown 
in Figure 6.1. Note that this type of instability is observed in crystal growth 1,2 
and in laminar flame propagation. 3,4 

We now shall study the nonlinear evolution of unstable modes. A general 
approach to this problem is based on the numerical solution of heat conduc¬ 
tion equation for the condensed phase with appropriate initial and bound¬ 
ary conditions. 5,6 A numerical solution of the three-dimensional heat con¬ 
duction problem with nonlinear boundary conditions at a moving boundary, 
whose shape and position were not known beforehand, requires a consider¬ 
able amount of computer time. The numerical calculations 5,6 are limited to 
studying the two-dimensional problem in which the temperature field is inde¬ 
pendent of the ^-coordinate. This approach reduces considerably the volume 
of numerical calculations while preserving some important features of the non¬ 
linear development of the corrugation instability. The equations to be solved 
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Figure 6.1: Dispersion curve and wavenumbers of the initial perturbations. 


are those considered in Chapter 5 (see Equations (5.1) through (5.3)). For 
convenience, we write them in the dimensionless form as: 

dO/dr = d 2 0/d £ 2 -f d 2 0/dr} 2 + exp(Z - Q 

dO/d C - dO/dr) • dZ/drj = XZ(1 - A K) 

Z = v [l — (dZ/drj) 2 ] 1 ^ 2 ex p [—(1 — A K)p/0] 

where 

O — finTjAI £ = fiz r) = fiy \ — finL v /cAI 

T = t/xn 2 v — Vo/fix A = an/pL v K = l/fiR 

R is the radius of curvature of the perturbed evaporation front and the func¬ 
tion £ = Z(rj, t) defines the position of the phase boundary. The dimensionless 
curvature of the phase boundary, K is given by: 

K = (d 2 Z/dr?) [1 + (dZ/dr)) 2 ] ~ 1/2 (6.3) 

The problem (6.1) through (6.3) is best investigated numerically in a labora¬ 
tory framework, rather than in the framework associated with the vaporization 
wave front. We consider periodic perturbations of the temperature and phase 
boundary shape: 

0(v,C,t) = 0(r) + 2n/k,C,T) Z(t},t) = Z(r) + 2n/k,T) 

where k is the perturbation wavenumber. The initial conditions to Equations 
(6.1) through (6.3) are chosen as: 

0 ( 77 , C,0) = 0 and Z(rj,0) = Zo{rj) (6-4) 

The problem formulated above was solved by numerical methods. 5-7 The de¬ 
tails of the numerical procedure can be found in Gorberg. 7 Here we will con¬ 
sider the results. First, we examine the situation corresponding to a slightly 


( 6 . 1 ) 

( 6 . 2 ) 
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supercritical laser intensity: I > I* and I — J* J*. For this case, a small 
interval of wavenumbers, &i 7 < k < k 2l , £27 — fci 7 fci 7 , exists in which the 
linear growth rate 7 (fc) is positive. The computations 5 show that when the 
initial conditions have the form of “white noise” i.e., the function Zo(rj) con¬ 
tains a large number of Fourier harmonics of approximately equal amplitudes, 
then the nonlinear evolution of instability is universal, regardless of the spe¬ 
cific initial shape of the phase boundary. At the initial stage, when the mode 
amplitudes are small, the growth of each mode is determined by the linear dis¬ 
persion equation (5.11); with increasing time, the amplitudes of the unstable 
modes (with k\ 7 < k < k 27 ) increase exponentially, and the amplitudes of the 
modes lying outside the instability interval (fci 7 ,fc 2 7 ) decrease exponentially. 
As the amplitudes of the unstable modes increase, the resonant excitation of 
higher harmonics takes place. The amplitudes of these harmonics also increase 
with time, despite the fact that the wavenumbers of the harmonics lie in the 
region of linear stability. Since every mode modulates the temperature field 
with its own spatial period, the interaction between different modes (not as¬ 
sociated with resonant conditions ki — nkj where n is an integer) leads to the 
suppression of all modes, except the mode which has the greatest amplitude 
when nonlinear interaction of modes begins and its resonant satellites. 

For the initial perturbations of the type of “white noise” this is a mode 
with maximum linear growth rate (see Figure 6 . 1 ), since it grows faster than 
any other during the initial linear stage. Note that for low supercriticality, all 
the unstable modes are concentrated in a small vicinity of k m . For this reason, 
two unstable modes cannot satisfy the resonant condition ki = nkj , and one 
can expect that only two modes with k « k m and k « 2 k m will exist at r 1 . 
An example of computations illustrating the evolution of perturbations in a 
slightly supercritical case is presented in Figures 6.1 through 6.3. Figure 6.1 
shows the linear growth rate 7 (k) with the wavenumbers fc n , which satisfy the 
condition of periodicity. The initial condition for numerical computations was 
chosen as a superposition of modes with k = k n : 

N 

Zo(v) = 52 •‘‘M 0 ) sin ( 6 - 5 ) 

n —1 

where k n = nk\, N = 10, and A n (0) =0.0015 = constant. In Figure 6.2 the 
nonlinear evolution of a vaporization front is shown. The front shape may be 
represented in the form of the Fourier series: 

Z{n, t) = Ao(t) + 52 A n(r) sin [fc n + i/>„(r)} 

n 

where the amplitudes and phases of the Fourier components vary with time. 
The amplitudes A n are shown in Figure 6.3 as functions of time. We see from 
the computations that after some transient process a steady state nonplane 
vaporization wave moving at constant velocity is formed. The fundamental 
spatial period of the wave is equal to 27r/A;4, where k± is the wavenumber of 
the mode having the largest value of the linear growth rate (see Figure 6 . 1 ). 
From Figures 6.2 and 6.3, one can see that the shape of the front is not 
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Figure 6.2: The nonlinear evolution of a vaporization front. Initial pertur¬ 
bations of the type of “white noise”. Perturbation wavenumbers k n = nki. 
The values of r are indicated on the left. The distances between successive 
front positions are reduced by a factor of 10. The scale of the front structure 
remains unchanged. 



Figure 6.3: Fourier components of the structure shown in Figure 6.2 as func¬ 
tions of time. 
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sinusoidal. The resonant harmonic having the wavenumber k$ = 2k 4 has a 
large amplitude—despite the fact that this mode is located inside the linear 
stability region. 

Consider now the case of an almost monochromatic initial condition. As¬ 
sume that the initial condition (6.5) besides the “white noise” modes contains 
a mode with k = kj, where kj is notably different from k m . The amplitude 
Aj of this mode is supposed to be much larger than the amplitudes of the 
“white noise” components A n . The computations show that, if at the onset of 
nonlinear interaction between the modes the amplitude Aj is greater than the 
amplitude A r of the harmonic whose wavenumber k r is closest to k m , which 
corresponds to the maximum linear growth rate, the interaction of these two 
modes leads to the formation of a stationary structure with fundamental pe¬ 
riod 27T /kj. If, at the nonlinear interaction stage, the amplitude Aj is less 
than the amplitude A r of the mode whose linear growth rate is higher, the 
initial condition is “forgotten” and a structure is formed with a spatial period 
2ir/k r . This is illustrated in Figures 6.4 and 6.5, which show the evolution of 
a vaporization front and the time dependence of the Fourier amplitudes A n 
for the initial conditions (6.5) with A n = 0.0015 (n ^ 3), A 3 = 0.045. 

We see that the structure with the wavenumber k% is suppressed at a 
moment r « 50, despite the fact that its initial amplitude was fairly large and 
that its growth rate 7 (^ 3 ) is positive. After a certain transient (50 < r < 70), 
a stationary wave having a period 2n/k 4 is formed. We see, thus, that the 
structures having a fundamental period 2n/k, whose wavenumbers k lie within 
the instability interval (Aq 7 , & 2 7 ) but close to its boundaries, are unstable and 
disintegrate. To understand this fact, a closed phenomenological equation 
was proposed 5,8 that describes the shape of the phase boundary in slightly 
supercritical regime. It was shown that the equation 

dZ/dt + LZ + ( dZ/drj) 2 = 0 , ( 6 . 6 ) 

where LZ — d 4 Z/drj 4 -f 2 ad 2 Z/drj 2 -f Z, a = k^ = (7 m -f l ) 1 / 2 « 1 + 
7 m /2 and 7 m = 7 (k m ) is a maximum value of 7 (k) in the instability interval 
(^ 17 ^ 27 ) (see Figure 6 . 1 ), accurately depicts (in two-dimensional case) the 
results of numerical calculations, allowing a simple investigation of the general 
properties of unstable evaporation. 

Equation ( 6 . 6 ) has stationary solutions in the form: 

Z(rj) = Asinkrj A 2 = a 2 — 1 — (fc 2 — a ) 2 (6.7) 

where k lies in the interval of instability of the trivial solution Z(t),t) = 0 : 
[a — (a 2 -!) 1 / 2 ] 1 / 2 < k < [a-f (a 2 -!) 1 / 2 ] 1 / 2 . Equation ( 6 . 6 ) was employed for 
the analysis of stability of steady-state solutions (6.7). The stability criterion 
can be written in the form: 


A 2 (k) > 2A 2 m /3 ( 6 . 8 ) 

where A m = 7 m = (a 2 — l ) 1 / 2 is the maximum value of the amplitude of 
stationary solution (6.7). For a given supercriticality, from Equation ( 6 . 8 ), we 
see that only those stationary structures whose amplitudes are not too small 
are stable. Those structures whose k lie near the boundaries of the plane front 
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Figure 6.4: The nonlinear evolution of a vaporization front. The initial am¬ 
plitude of perturbation with k = k-j exceeds the level of “white noise” by a 
factor of 30. 



Figure 6.5: Fourier components of the structure shown in Figure 6.4. 
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Figure 6 . 6 : Isotherms of the temperature field in the stationary evapora¬ 
tion wave shown in Figure 6.2. Curve 1—0 = O.710 m , 2—6 = O.780 m , 

3—0 = O.840 m , 4 —6 = O.910 m , 5 —6 = O.970 m . 



instability range are absolutely unstable. Note that an analysis of stability of 
plane Couette flow leads to a similar criterion . 9 Using condition ( 6 . 8 ), we can 
conclude that the stationary solution, with fundamental spatial period 2 Tt/kz, 
is unstable and should be destroyed at arbitrary initial amplitude A 3 . 

We now shall discuss the characteristic properties of the stationary struc¬ 
ture. The temperature distribution in a nonplane steady-state evaporation 
wave with fundamental period 27r/A;4 (see Figure 6.1) is shown in Figure 6.6. 
The phase boundary is not isothermal, and those portions of the boundary 
displaced to the condensed phase have higher temperatures. Because the evap¬ 
oration wave is stationary and its propagation velocity is constant, the effective 
vaporization energy L e // in the “hot” portions of the front must be higher 
than in the “cold”. The curvature of the isotherms decreases quickly with dis¬ 
tance from the phase boundary. The isotherm with temperature 0 = O.70 m , 
where 0 m is the maximum temperature, is almost planar and its position is 
the same as in the case of the plane vaporization front corresponding to the 
same laser intensity. 


6.2 Strongly supercritical regimes 
of laser evaporation. 

“Spin” sublimation of solids. 

Previously we have discussed the case of low supercriticality when the interval 
of unstable modes is narrow and the second harmonic of the perturbed funda¬ 
mental mode lies in the region of linear stability. In this case, the evolution of 
the perturbations results in the formation of stationary structures with non¬ 
plane, nonisothermal phase boundary. As the supercriticality increases, the 
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wavenumber of the second harmonic approaches the upper boundary & 2 7 of the 
instability range and the situation may change qualitatively. If the intensity / 
of the laser radiation greatly exceeds the instability threshold, the range where 
the linear growth rate is positive can become so wide that multiple harmonics 
of the monochromatic initial perturbation will lie in it. 

This situation occurs when the laser intensity I is greater than I* deter¬ 
mined from the equation 2 &i 7 (J*) = k 2 ~ i {I *). Calculations show that in this 
case an important role is played by the generation of harmonics, which leads to 
nonstationary regimes of evaporation front propagation. Actually, the steady- 
state propagation can be violated even earlier, at I somewhat smaller than 
/*, since the coupling between the modes can compensate for the small damp¬ 
ing of the second harmonic. We refer to Figure 6 . 1 . Suppose that the initial 
perturbation has the wavenumber fc 3 , lying in the instability range (&i 7 , £ 27 )- 
The second harmonics of this perturbation has wavenumber fcg = 2k 3 located 
in the region of linear stability near the upper boundary £27 of the range of 
unstable modes. The damping of the second harmonic in this case is small 
and may be less than the amplification resulting from mode interaction if the 
amplitude of the fundamental mode k = k 3 is fairly large. One might expect 
a stationary regime of evaporation not to be reached at this level of supercrit¬ 
icality. A numerical computation 6 shows that, after a certain transient, an 
asymptotic traveling-wave mode is established. We recall that in Anisimov et 
al . 6 a two-dimensional temperature field was considered. A physical example 
related to this calculation is the vaporization of a thin-walled cylinder irradi¬ 
ated along its axis. Here, the radial component of the temperature gradient 
may be neglected, so that the temperature 6 = 0(rj, C, r) where the 77 coor¬ 
dinate is measured along the circumference of the base of the cylinder and 
the £ axis is parallel to the cylinder axis. The ( 77 , Q plane may be considered 
as the convolution of the side of cylinder. The radius, #, of the cylinder is 
determined by the condition of periodicity: 

6 (rj -f 2ttR, £,t) = 0(77, £, r) or R=l/k\. 

With this constraint in mind, we now can interpret the results of computations 6 
as the formation of the structure which rotates around the cylinder axis and 
propagates with constant velocity along this axis. The points on the front 
having the same temperature form helical lines winding along the side of the 
cylinder. Similar spinning structures were observed in combustion waves . 10 In 
Figure 6.7 steady-state propagation of the above described “spin” evaporation 
wave is shown. The steady-state regime is reached after a certain transient. 
The duration of the transient and the direction of wave motion along the 
77 -axis are determined by the specific initial conditions. 

In the Fourier expansion of the phase boundary shape, the modes with 
wavenumbers k 3 and k& = 2k 3 are present. The interaction between these 
modes leads to a periodic variation of their amplitudes. The law of motion of 
the vaporization front is described by the following equation: 

Z(rj, r) = (3 t + A 3 sin(fc 3 77 - (3t) + A 6 sm(k 6 r) - (3t + \j ) 0 ) 
where /?, A 3 , Aq , and t/>o are constants. 
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Figure 6.7: Steady-state propagation of the “spin” 
vaporization wave. The displacement of the front 
is reduced by a factor of 6 . The time intervals be¬ 
tween subsequent locations of the front At = 11.4, 
l — 27 vR = 2n/k3. 



When the laser intensity I reaches /*, the second harmonics wavenumber 
falls inside the range of linear instability. Numerical computations show that 
when this happens, a complex, nonstationary propagation mode arises. An 
example of this mode is presented in Figure 6 . 8 . Unlike the case of “spin” 
evaporation, a Fourier series expansion of the function Z(r),r) in this case 
shows a nonperiodic time behavior of the front shape. However, to answer 
the question of whether this mode is complex periodic, quasi-periodic, or true 
stochastic, however, one needs more information concerning the lasting behav¬ 
ior of Fourier amplitudes. We note that this question seems to be academic, 
because in experiments it is practically impossible to establish decisively which 
of the complex, nonstationary vaporization modes occurs. 


6.3 Nonlinear stages of instability 

of evaporation waves in dielectrics 

As was shown in Chapter 5, there are two different modes of instability of 
plane evaporation waves in dielectrics. The first is similar to the corrugation 
instability of evaporation waves in highly absorbing solids. For this mode the 
growth of perturbations is aperiodic and the growth rate reaches its maximum 
at finite k. The second one is the instability mode similar to thermal explosion. 
Maximum growth rate is reached for this mode at k = 0 . When J£o 1, 
there is no solution of the thermal-explosion type. If k £o 1, the growth rate 
is given by Equation (5.38): 

7 (fc) = - JV S + v.(k 2 - JV s /x) 1/2 (6.9) 

from which we can see that the perturbations may grow in an oscillatory 
manner if k is small. A more detailed analysis shows that oscillatory solutions 
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Figure 6.8: The shape of the phase boundary for the “stochastic” vaporization 
mode. The distances between subsequent locations of the front are reduced 
by a factor of 20. The values of r are indicated on the left. 
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exist for arbitrary J. However, the range of wavenumbers corresponding to 
these solutions decreases rapidly as J increases. 

The region where oscillatory solutions exist is shown cross-hatched in Fig¬ 
ure 5.4. It is seen that in a wide range of laser intensities, the oscillatory 
solutions are possible if £o is sufficiently small. One may expect that the 
nonlinear evolution of some modes from this range will lead to the formation 
of plane oscillatory evaporation waves. This assumption is confirmed by nu¬ 
merical calculations. 11,12 The calculations 11 were performed for a simplified 
model in which the temperature dependence of thermal conductivity was not 
considered, and the light absorption was described by the simplest Bouguer 
law. The results of calculations 11 show that three types of plane vaporiza¬ 
tion waves exist in nonlinearly absorbing media: (1) stationary waves stable 
to perturbations with k = 0; (2) oscillatory waves; and (3) waves in which 
vaporization ceases because of bleaching of material. In reference 12 , the elec¬ 
tronic heat conduction was taken into account and the light absorption was 
calculated by solving Maxwell’s equations for the electromagnetic field of a 
light wave. We now discuss these calculations more in detail. In this work, as 
in Chapters 1 and 2, the energy transfer in condensed phase is described by 
the heat conduction equation 

cpdT/dt = d/dz[(n o -f n e )dT/dz ] -f cpVdT/dz + Q (6.10) 

where ko and n e are the phonon and electron components of thermal con¬ 
ductivity. The heat source, Q, in Equation (6.10) is given by formula (1.24). 
To determine the electric field the wave equation (1.25) is solved numerically. 
The complex dielectric constant of ionized dielectric is written in the form 
(see, Poyurovskaya et al. 13 ): 

e = e 0 — Ama/u 6 0 = tIq — 47T a jv and a = gqv 2 /(v 2 + lUq) 

where a is the electrical conductivity at laser frequency u>o; is the static 
conductivity; v is the effective collision frequency; and no is the refraction 
index of the cold dielectric. The thermal and static electrical conductivities 
are assumed to be proportional to the concentration of free electrons, that is 
strongly temperature dependent: 

k c (T) = ttgoexp (-Eg/T) and cr 0 (T) = a 0 oexp(-E g /T) (6.11) 

The usual boundary conditions (1.18) are imposed at the vaporization front 
z = Z(t). 

Equation (6.10) has a steady-state solution T s (z). This solution was 
found 12 by numerically solving Equation (6.10) (with dT/dt = 0) together 
with the wave equation (1.25). Then, a small perturbation was introduced and 
its evolution was studied by solving numerically the complete time-dependent 
equation .(6.10). The computations show that in a certain range of laser in¬ 
tensities, the stationary solution T s (z) is stable. Initial perturbations of the 
temperature field decrease with time, so that the temperature and the front 
velocity tend to their stationary values: T(z,t) —> T s (z), and V(t) —> V s . 
Outside the stability region small one-dimensional perturbations grow with 
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Figure 6.9: The self-excited oscillation mode of dielectric vaporization. Depen¬ 
dence of the front velocity v and the reflection and transmission coefficients, 
R and T, on time. 


time. The evolution of the perturbations and the character of the asymptotic 
propagation mode depend on the laser intensity and parameters of the mate¬ 
rial irradiated. Three asymptotic regimes were observed in computations: 12 
(1) stable nonlinear oscillations with periodic variation of the front velocity 
(Figure 6.9); (2) low-temperature regime with complete bleaching of a speci¬ 
men (Figure 6.10); and (3) separation of the absorbing layer from the phase 
boundary and formation of ionization wave not connected with phase transi¬ 
tion (Figure 6.11) Note, that the third regime is possible only when light is 
incident from the side of cold dielectric. 

As is seen in Figure 6.9, the oscillations of vaporization front velocity 
are accompanied by the oscillations of reflection, R , and transmission, T, of 
laser radiation. The period of oscillations is, in order of magnitude, equal to 
x{T S m)/Vs > where x(T S m) is the thermal diffusivity at the maximum steady- 
state temperature, T 5m . The radiation is incident on the specimen from the 
side of vaporization front. A similar regime also exists when the radiation 
is incident from the side of the cold dielectric. The parameters of the tar¬ 
get material in the calculation presented in Figure 6.9 are close to those of 
fused quartz. For these parameters of the target and the same laser intensity 
there exists, formally, a stationary solution of Equation (6.10) with tempera¬ 
ture profile T s (z) and front velocity V a , which is unstable. In the stationary 
evaporation mode, the light-absorbing layer (i.e., the neighborhood of the 
point where the temperature reaches its maximum) and the phase boundary 
move at the same velocity, V s . In the oscillatory mode, the distance between 
these regions changes periodically with time. Over some regions of parame- 
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Figure 6.10: Unstable modes of dielectric vaporization. Irradiation from the 
side of the vaporization front; the front moves more rapidly than the absorbing 
layer. 



Figure 6.11: Unstable modes of dielectric vaporization. Irradiation from the 
side of cold dielectric; the absorbing layer moves faster than the vaporization 
front. 
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ters, the synchronism in propagation of the absorbing layer and vaporization 
front may break down. If the vaporization front moves faster, it overtakes 
and ” swallows” the absorbing layer. The optical thickness of the absorbing 
layer rapidly decreases and the specimen becomes transparent. An example 
of one-dimensional instability leading to the bleaching of the dielectric tar¬ 
get is shown in Figure 6.10 (the laser beam comes in from the side of the 
vaporization front). 

A different mode of instability may arise if the absorbing layer moves 
faster than the vaporization front. In this case, the absorbing layer breaks 
away from the phase boundary, and the heat flux to the vaporization front 
decreases. This leads to the deceleration of the phase boundary and formation 
of an absorption wave similar to the optical breakdown waves in gases. 13,14 
According Anisimov et al., 12 this instability occurs when the laser radiation is 
focused on the exit surface of dielectric sample (irradiation from the side of the 
cold dielectric). An analysis of the numerical results shows that the evolution 
of the perturbations of a plane stationary evaporation wave is determined by 
two parameters: (1) the ratio P = U/E g of characteristic activation energies 
for evaporation ( U ) and for high-frequency conductivity ( E g ), and (2) the 
dimensionless laser intensity Q = Io/pLoVo. If radiation is incident from the 
vaporization front side, and the parameter P does not exceed some critical 
value Pi « 1, then as Q increases at constant P, the following changes occur: 
the stationary evaporation changes into a self-sustained oscillation mode and 
then to specimen bleaching due to the absorption layer being vaporized. A 
further increase in the laser intensity leads to a change in propagation modes 
in the reverse order. If P > Pi, a stationary evaporation is stable with 
respect to perturbations with k = 0. When light is incident from the side of 
the cold dielectric, and P < P 2 « 1, the sequence of modes as Q increases 
is: separation of the absorbing layer from the phase boundary, bleaching of 
the target due to vaporization of the absorbing layer, oscillatory mode, and a 
stationary evaporation wave. 
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Chapter 7 

INSTABILITIES IN THE 
LASER INTERACTION 
WITH LIQUIDS 


7.1 Thermocapillary instability of a liquid 
heated by laser radiation 

In this section we will examine the instability of a plane homogeneous liquid 
film which absorbs laser radiation. The instability arises from the temper¬ 
ature dependence of surface tension <t(T). The capillary forces are small in 
magnitude; therefore, to observe the liquid motion induced by them it is nec¬ 
essary that other forces—in particular, reactive forces which arise during the 
evaporation of liquid—do not act on the layer. Therefore, we would expect 
an instability of this type to be observed at relatively low laser radiation in¬ 
tensities in materials whose vapor pressure at melting point is low. 

To understand qualitatively how instability arises, we shall examine a liq¬ 
uid film of thickness, H , situated horizontally in a gravitational field and 
heated from above by laser radiation. We shall consider the radiation to be 
absorbed in infinitesimally thin surface layer because, as will be seen, long¬ 
wave perturbations have the greatest growth rate. For simplicity, we shall 
assume that the temperature gradient close to the surface is constant in the 
steady state: 

dT^Idz — AIjn— constant (7.1) 

Here, I is the laser radiation intensity, and A is the fraction of laser intensity 
absorbed by the liquid. 

It is well known 1 that there are two types of oscillations which can be 
excited in liquid film under examination: gravitational-capillary and ther¬ 
mocapillary. Excitation of thermocapillary oscillations is associated with the 
temperature dependence of surface tension. Because the surface tension de¬ 
creases with temperature, 2 forces will occur on a nonuniformly heated liquid 
surface that are directed away from the heated portion toward the cold por¬ 
tions and will cause motion in the layer of liquid adjacent to surface. If the 
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Figure 7.1: Flow structure caused by thermocapillary force. The dashed line 
represents the temperature distribution; the solid lines represent the flow. 


surface temperature profile is periodic, a liquid flow will be set up as shown 
in Figure 7.1. It can be seen from the flow pattern that liquid in the surface 
layer will flow from the heated portions to the cold portions and that liquid at 
a lower temperature will rise from below to replace it. Convective transport 
and thermal conduction will equalize the temperature and, due to the inertia 
of the liquid, the equalization process may have an oscillatory nature. 

It follows from the arguments presented here that, should the liquid tem¬ 
perature not fall, but increase from below, the thermocapillary effect may lead 
to an aperiodic instability. Such an instability was studied by Pearson, 3 and 
will not be considered here. 

If the gradient of temperature of the liquid has the sign defined by Equation 
(7.1), the situation becomes more complex. The interaction of two branches of 
spectrum of the natural modes—gravitational-capillary and thermocapillary 
oscillations—may lead to a parametric amplification of the fluctuations, i.e., 
to the occurrence of instability. For a description of this phenomenon, we 
shall introduce a scalar potential p and a vector potential A of the velocity, 
satisfying the equations: 

v = V<p -f rot A div A = 0 

where v(r,£) is the velocity of liquid. Introducing a temperature perturba¬ 
tion Ti(r, t) = T(r, t) — To(z), \T\\ To, and linearizing the Navier-Stokes 
equations and the heat conduction equation for small values of v and T \, we 
obtain the set of equations: 



dA 


dt 

dT\ 

dt 


— v AA 


+ v s 


dT 


o 


dz 


0 

0 

xAT x 


(7.2) 


where v = rj/p is the kinematic viscosity and x = ^/cp is the thermal diffu- 
sivity of the liquid. The boundary conditions at the lower boundary of the 
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film, x = —H have the form: 


T\ = 0 and v = 0 



On the upper boundary, the boundary conditions have to be imposed on the 
perturbed liquid surface z = Z(x,y,t). In such cases, it is convenient to 
project the boundary conditions onto the plane z = 0 by expanding all the 
functions occurring in the power series in Z. Then, in the linear approximation 
with respect to Ti, v and Z, the boundary conditions to the Navier-Stokes 
equations (with the conditions that all three force components acting on unit 
area of the surface of a liquid must be equal to zero) have the form : 1 


d 2 ip 

p -dW- a 



d 2 v z 

+ P9V ‘ + 2v Mi 



da d f dT\ dTo 
dT dx V dt Vz dz 



dvx 

dz 


+ 


dv z 

dx 


da d 7 dT\ dTo \ d / dv y dv z \ 

dT dy \ dt + Vz dz ) ^ dt \ dz + dy ) 



In general, the boundary conditions for the heat conduction equation, that 
ensures the continuity of normal component of energy flux through the phase 
boundary, for this problem may be written in form: 


— ndT/dn + AI = 0 when z = Z(x,y,t) 


where d/dn is the normal derivative. Projecting this onto the plane z = 0, 
differentiating with respect to time, and taking into account the fact, that 
dZ/dt — v z , we obtain 


d 2 T x , (d 2 Tp\ 
dtdz Vz V dz 2 ) 


when z = 0 


However, according to our assumption (see Equation (7.1)), d 2 To/dz 2 = 0 , 
from which it follows that 


dT\/dz = 0 when z = 0 


(7.5) 


An arbitrary smooth perturbation to the shape of the liquid surface and 
its temperature may be expanded in the Fourier integral. Therefore, when 
using the principle of superposition in the linear approximation, to determine 
the evolution of an arbitrary perturbation, it is sufficient to investigate the 
evolution of the monochromatic perturbation of the form: 

ip = <p(z) exp (ikx -f 7 1 ) A y = A(z) exp (ikx + 7 1 ) 

(7.6) 

T\ = Ti(z) exp(zA;a: + 7 1 ) 
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Substituting (7.6) into (7.4), we obtain: 



0 

0 




-f ikA(z) 


q = k 2 +i 

X 


(7.7) 


(7.8) 


From the condition that Equation (7.7) should have a solution satisfying the 
boundary conditions (7.5) and (7.6), the dispersion equation 7 (k) follows for 
the perturbations considered. However, the investigation of Equation (7.6) is 
so difficult for a general case that we will confine our analysis to the case of 
a “deep liquid” kH 1 . A more general study is performed in the original 
work . 4 The dispersion equation, when kH 1 , becomes 


- N('y,k)M('y,k)] + 2 N('r,k) = 1 + (^ (7.9) 

where 

u>$ = gk + (a/p)k 3 M{~i, s) = \h ~ l/h(h + s) 

k) = (u>§ + 2 -yvkli) / (u>$ + j 2 + 27 vk 2 ) 

and l\ and I 2 are defined by Equation (7.8). From Equation (7.9) in the 
limiting case when (x^ 2 **^ 2 ) (I 7 I ,^o) we easily obtain 


where 


( 7 2 + w o)( 7 2 + c 2 k 2 ) = 0 



da dTo 
dT dz 


p -l(l + P l/ 2 )-l 


(7.10) 


and P = vjx is the Prandtl number. From Equation (7.10), it can be seen that 
the dispersion equation has two independent oscillatory branches in the long¬ 
wave limit: gravitational-capillary, with 7 ^ = iuJo(k) and thermocapillary 

with = ick. Assuming further that 7 i(fc) = iw${k) -F £i(fc), 72 {k) = 
ick + S 2 {k), the relation between modes leading to parametric instability may 
be obtained. Instability will begin if Re<5i ,2 > 0 - The calculations show 
that the sign of Re 72 is the same as that of the difference ck — As 
is clear from Equation (7.1), the magnitude of c is proportional to (A/) 1 / 2 , 
and u;q is independent of the radiation intensity. This means that a certain 
critical level of intensity I * exists, above which Re £2 becomes positive, i.e., 
thermocapillary waves are excited. From Figure 7.2 it is clear that the wave 
numbers of the excited modes lie in the finite range of k\(I) <k <k 2 {!)■ 
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1 



Figure 7.2: Dispersion curves Im 7 (fc). 1 —Gravitational-capillary waves; 

2—Thermocapillary waves, I> I*; 3 —1 = 1*; and 4—/</*. 

The behavior of Re£i is somewhat more complex: it is negative at u>o < ck , 
but when u )o > ck, for every k there is a certain minimum I m (k) above which 
Re£i > 0 . Calculations show that: 

mm[/ TO (A:)] = I m (k c ) = 2 5/2 t x) l,2 {k c A{da / d,T))~ l 

where 

k c = (pg/5a) 1/2 and u c = LJ 0 {k c ) 

We now shall consider the case of resonant interaction of natural modes. 
This case corresponds to those values of k which are close to the solutions 
of the equation u>o(A;) = kc. Near resonance the real part of 7 increases and 
ceases to be a small correction to the imaginary part. Calculations show that, 
for this case, 


Re 7 ~ (o; 0 /2)(xA: 2 /^o) 1/4 (1 -h P 1/2 ) 1/2 sin( 7 r/ 8 ) 

Whereas, far from resonance, it is possible to distinguish potential (gravita¬ 
tional capillary) and vortex (thermocapillary) modes, such a distinction be¬ 
comes impossible in the resonance regions. In the developing complex motion, 
both components—potential and vortex—have the same order of magnitude. 

The analysis, performed by Levchenko and Chernyakov 4 for a liquid of 
finite depth, shows that the qualitative picture of the occurrence of instability 
described above is preserved completely in the general case. The quantitative 
difference is in the increase in the threshold intensity /*, which is determined 
by the additional dissipation of energy near the bottom boundary of the liquid 
film. 

We will briefly discuss the nonlinear stage of the instability development. 
From the discussion above, it is clear that the liquid remains motionless, and 
the temperature profile To{z) is linear for radiation intensities I < I min = 
Im{kc) m For a small excursion above the threshold: I > I mini I ^min 
I min , motion whose characteristic scale is A c = 2ir/k c occurs in the liquid 
film. It is natural to assume that the growth of the unstable mode in the 
weakly supercritical case will stabilize due to nonlinear effects; this leads to 
the formation of a new stationary state with broken symmetry. The spatial 
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structure and the temperature and velocity fields in this state may be different 
and, for given k c , there may be rolls, as well as triangular, rectangular, or 
hexagonal cells. 

Finally, it is possible, generally speaking, to have a case where stationary 
structures are not formed at all, and the motion of the liquid has a complex, 
chaotic, turbulent form. Of all the possibilities enumerated by Levchenko 
and Chernyakov 4 only the simplest two-dimensional structure was studied 
in the form of rolls. If the supercriticality is small, we can assume that an 
almost monochromatic gravitational-capillary wave occurs whose amplitude 
Zk satisfies the equation: 



2 

However, 7 itself depends on \Zk\ in the nonlinear case and may be presented 
in the form of the expansion: 

7 = Re ^1 4- P | Zk | 2 + ... 

To calculate the coefficients of this expansion, standard methods of bifurcation 
analysis 5 are used: the unknown functions are expanded in a formal series 
in powers of e = kZ fc, where Zk is the surface displacement amplitude. In 
performing such calculations, it is convenient to convert to new independent 
variables : 6 the potential <p and the stream function </>• The result obtained 4 
is that the perturbed mode is soft, and its amplitude is 

Zfc ~[/-/ m (*)] 1/2 

when I > I m i n . The stability of this structure was not investigated. 

Strictly speaking, it follows from the study of weakly supercritical struc¬ 
tures that the dependence of the surface reflection on the curvature of the sur¬ 
face must be considered (even for normal incidence of the light). The results 4 
were obtained neglecting this dependence. However, there is no reason to ex¬ 
pect that taking into account radiation diffraction by the perturbed surface 
will qualitatively change the results, especially the conclusion concerning the 
soft character of wave excitation. 

To conclude this section, we present numerical estimates, which show the 
orders of magnitude of the fundamental quantities characterizing the insta¬ 
bility being examined. For the melting of iron, the threshold value of the 
absorbed intensity is Af m i n ~ 10 3 W/cm 2 . For the intensity considerably 
higher than the threshold value, AI = 5 x 10 5 W/cm 2 , the resonance wave 
number equals k ~ 500 cm -1 , and the corresponding value of the instability 
growth rate equals Re 7 ~ 5 x 10 4 sec -1 . 


7.2 Corrugation instability of plane 
evaporation waves in liquids 

Effects of instability which destroy the simplest pattern of evaporation with 
a plane front are especially important if a liquid layer is formed on the solid 
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body surface. As was shown in Chapter 1, this situation is typical for evap¬ 
oration of metals. Experimental studies show that even with rather uniform 
laser intensity profile, a plane evaporation front turns out to be unstable. 7,8 
The distortion of a plane evaporation front can be attributed in many cases 
to the material inhomogeneity, spatial fluctuations of laser radiation intensity, 
and volume vaporization (boiling). 9,10 These effects can be eliminated to a 
large extent by a proper choice of experimental conditions. A more general 
cause of liquid-vapor phase boundary distortion is the intrinsic instability of 
vaporization front that moves in the direction of temperature gradient. 11,12 In 
Chapter 5 this instability was studied for a solid-vapor phase transition. In 
this chapter, we consider the interaction of a spatially uniform laser beam of 
constant intensity with an absorbing liquid layer, and investigate the stability 
of surface evaporation of a liquid under these conditions. 13,14 We assume that 
the gravity force (natural or artificial) is directed normal to the liquid surface. 
We neglect the flow of liquid along the surface of the solid due to the nonuni¬ 
formity of the reactive pressure, the boiling of liquid, and the temperature 
dependence of the parameters of liquid. We take into account, however, the 
finite thickness of the liquid layer. 14,15 

As in the case of sublimation, it is convenient to choose a coordinate frame 
in which a plane stationary evaporation front is at rest and lies in the plane z = 
0, the axis z being directed toward the condensed phase. In this framework, 
the condensed material moves along the z-axis with the stationary velocity 
v z = — V s . As in Chapters 1 and 5, the temperature field in the condensed 
material is described by the heat conduction equation: 

cp ( dT/dt + v • VT) = div (/cVT) + pi (7.11) 

with the boundary condition T —> 0 when z —► oo and the energy balance 
equation on the evaporation front at z = Z(x,y,t) 

L eff V nP = K(VT) n (7.12) 

The evaporation rate is related to the surface temperature: 

V n = Vo exp(—ML e/ f /k B T) (7.13) 

The effective heat of vaporization depends on the local curvature of the phase 
boundary as (see Equation (5.3)): 

L e /f = L v — (a/p)(l/R\ + I/R 2 ) 

where R\ and R 2 are the principal radii of curvature and a is the surface 
energy density. 

A number of interesting features of laser evaporation may be attributed to 
the diffraction of monochromatic laser radiation on the relief formed by the 
disturbances of the evaporated surface. This effect is important (see references 
16 and 17) if the perturbation wavelength is close to the laser radiation wave¬ 
length. Far from this resonant region, diffraction effects may be neglected, 
and the intensity of laser radiation may be described by: 

I = Io exp [/i(Z - z)\ 


(7.14) 
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The normal component of the evaporation front velocity is determined by the 
vaporization and the hydrodynamic motion of the liquid: 

Vn = Vn d" (7.15) 

where u n is the projection of the mass velocity of liquid onto the normal 
to the phase boundary. The heat conduction equation (7.11) and boundary 
conditions (7.12) and (7.13) are quite similar to those used in Chapter 5. In 
the presence of a melt layer on the surface of the evaporated sample, the 
boundary conditions on the evaporation front and in the sample depth should 
be supplemented by the boundary condition on the melting front. In the 
case of metals, the density and heat capacity of the liquid phase are close 
to those of the solid phase. In our analysis we shall disregard the jumps in 
thermophysical properties of the condensed phase at the melting front. We 
also shall disregard the latent heat of fusion, which is much smaller than the 
latent heat of vaporization. These simplifying assumptions, normal for laser 
evaporation studies, lead to some distortion of the stationary temperature 
profile; but, they actually do not affect the results of stability analysis. Thus, 
in our formulation of the stability problem, among all the differences in the 
physical properties of liquid and solid phases only one property, namely the 
fluidity of liquid phase, which is most important for instability development, 
is taken into account. 

The melt may be considered, with sufficient accuracy, as an ideal incom¬ 
pressible liquid. With the purpose of linear analysis, we can restrict ourselves 
to the consideration of two-dimensional perturbations dependent on the co¬ 
ordinate z and one of the coordinates in the plane of the unperturbed evap¬ 
oration front (for certainty, y). The melt flow may be described by a stream 
function $(y,z,£) that satisfies the equation 

dA$ 3$ <9A$ 9A£_ o 

dt dz dy + dy dz (7.16) 

u z = d$/dy u y = —d$/dz 

On the evaporation front, the change in pressure is attributed to the effect 
of surface tension, p\ = a{\/R\ -f l/# 2 ), and to the recoil pressure of the 
evaporated material, P 2 = bpc s V n , where c s is the sound speed in vapor near 
evaporation front and the factor b & 1.67. The factor b takes into account the 
hydrodynamic boundary conditions on the evaporation front (Anisimov 18 and 
Chapter 3). The conditions of mass and momentum flux continuity are the 
boundary conditions for the hydrodynamics equations on the melting front. 
The continuity condition for the tangential component of velocity follows di¬ 
rectly from these conditions. 11 Emphasize that the condition imposed on the 
tangential component of velocity is in no way attributed to viscosity and is 
retained when considering inviscid fluid. If one neglects the difference in den¬ 
sities of the solid and liquid, the condition of mass flux continuity is reduced 
to the continuity of normal velocity component. In the presence of melting 
front perturbations, the continuity conditions for both, normal and tangential, 
components of the velocity can not be satisfied in assumption of a potential 
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flow. Thus, the melt flow turns out to be vortical. The boundary conditions 
of velocity continuity on the melting surface are equivalent to the conditions 
of continuity and smoothness of a stream function. 

To study the stability of a plane stationary evaporation wave, one should 
linearize Equations (7.11) through (7.16) in terms of small perturbations of 
temperature profile, stream function, evaporation rate, and the shape of the 
evaporation front. By virtue of the superposition principle, it is sufficient 
to consider disturbances proportional to <f>k{z) exp(iky + 7 £), where k is the 
perturbation wave number and 7 is the instability growth rate. After the 
linearization, the system of ordinary differential equations for the functions 
<j>k(z) is obtained. Solving these equations (with allowance for the boundary 
conditions on the melting front), we arrive at the system of three equations 
for the amplitudes of perturbations. It is convenient to write these equations 
in a matrix form : 14,15 


Aip = 0 (7.17) 

where p is the column of the amplitudes of evaporation rate perturbations, 
deformation of the evaporation front, and the mass velocity perturbation, and 
the matrix A is 


A = 


1 -f mc 2 xw/V 


bc.k 


Z 

7 1 

— k(g -f ak 2 /p) K 


(7.18) 


where m = cM/fcg, e = kBT s (0)/ML v , and uj is the root of the equation 


U) 2 - Vu)/x ~ (k 2 - 7/x) = 0 


(7.19) 


satisfying the condition Rea; > 0; 

-7 /■, , ,m/ 12 / \ V(xu> - V)(xu - V - mew) 

z = (i + ™x“V)Vok /p\ - x(u-n)-v - 

__ (i k 2 V 2 - 7 2 ) cosh(fcff) _ 

kV[cosh(kH) — exp (—7 H/V)\ — 7 sinh(A:i7) 

_ KV( 1 + me)xpP(p) - (V + mexp)P(V/x) 

2x(XM - V) cosh (kH) 

p, 0{k,p) 0(~k,p) _ 2kVe(-y/V,p) 

W (kV + 7) (kV - 7) ( k 2 V 2 - 7 2 ) 

= exp (tH)F(p + t + u> -V/x) F(z) = [1 - exp (~zH)]/z 


H is the thickness of the unperturbed melt layer, and V = is the stationary 
velocity of evaporation front. In the limiting case of significant melt depth, 
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kH 1 , for not very quickly damping solutions, Re ( 7 ) > -kV, Re^/V 4 
k)H > 1, 


K w kV 4 7 and $ 
and, for H —> 0 


V(k — id — V/x ~ mep) 


X{k 4- v)(k + id — V/x 4- p 

K « 2F/A;fl r2 and $ « FiJ/ 6 * 

The conditions for the existence of a nontrivial solution of Equation (7.17), 


det(A) = 0 , 


(7.20) 


is the dispersion equation for the instability growth rate 7 (k). 

Before studying the dispersion equation, it is useful to consider the struc¬ 
ture of the matrix (7.18). The upper left minor of the determinant (7.20) 
corresponds to the problem of pure evaporation (without the melt flow) and 
describes the sublimation of solids considered in Anisimov et al . 11 and in 
Chapter 5. Setting this minor equal to zero we obtain the equation: 

[l 4 me 2 (xid/V)] <y + Z = 0 (7.21) 

which, in combination with Equation (7.19) forms a set of equations equivalent 
to that obtained in Chapter 5. 

The bottom right minor of the determinant (7.20) describes the evolution 
of surface waves in a melt layer at a constant rate of melt evaporation. Sta¬ 
tionary liquid flow from the melting front to the vaporization front leads to 
significant differences of this case from ordinary gravitational-capillary waves. 
Having equated this minor to zero, we obtain the dispersion equation: 

7 K 4- (g 4 crk 2 /p)k = 0 (7.22) 


In the case of large melt depth it is reduced to a quadratic equation: 

7 2 4 kV 7 4 (g 4 ok 2 /p)k = 0 

whose solution is 

7 = —kV/2 4 \/[(kV/ 2) 2 — (g 4 crk 2 /p)k] 

The mass flux to the evaporation front from the melt depth results in the at¬ 
tenuation of gravitational-capillary waves. When g + crk 2 /p < 0 , the aperiodic 
growth of perturbations occurs. When kV 2 /4 > g 4 crk 2 /p > 0, a new mode 
of the evolution of perturbations—aperiodic attenuation—appears. With the 
decrease of kH , the oscillatory nonpotential modes become important. A typ¬ 
ical form of the dispersion relation for this case is shown in Figure 7.3. The 
above specific features of the dispersion dependence (7.22) remain also with 
the complete dispersion equation (7.20). 

It is convenient to write the dispersion equation (7.20) in the form: 

(7 + W)K + (-VW + g e /f + 9 + ok 2 /p) k = 0 (7.23) 

where 

g eff = (bc s + V)W and W = (Z + 7 $)/(l + me 2 X u/V - 3>) (7.24) 
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Figure 7.3: The dependence of the instability growth rate 7 on the pertur¬ 
bation wave number k and effective gravity acceleration (g ■f jk 2 ) for liquid 
film evaporation with constant evaporation rate. Solid lines—Re 7 ; dashed 
lines— Im 7 . 


It follows from (7.24) that, at k > V/x, W turns out to be close to — 70 , 
where 70 is the solution of Equation (7.21) corresponding to the case of pure 
sublimation: 

7o = ~Z/( 1 + me 2 xu>/V). 

The analysis of Equation (7.23) shows that W depends mainly on k. With 
the variation of other parameters ( g , H) there may occur strong changes in 7 
and relative variation of W is usually small. We may assume then approxi¬ 
mately that W « — 70 . Since the sound velocity c s is much larger than the 
vaporization front velocity V, the main effect of thermal perturbations at large 
k and H is described by g e ff- The sublimation growth rate thus describes 
the effective mass force affecting the melt. When there is an instability in 
the sublimation problem (70 > 0 ) then aperiodic instability may occur in the 
complete problem (7.23). The characteristic value of growth rate is greater 
than 70 : 

7 « ( 6 c s fc 7 o) 1/2 

The instability region is determined by the inequality: 

-geff « bc s 70 > g + ak 2 /p 

instead of 70 > 0 , and becomes, therefore, more narrow than in the case of 
sublimation. This mode of instability is suppressed completely at any radia¬ 
tion intensity when the (artificial) gravity acceleration is rather high, 

g > g c ~ max ((V 4- bc^kyo ~ °k 2 /p) 
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Figure 7.4: Dispersion curves 7 (k) at high radiation intensity J 0 = 5 x 10 7 
W/cm 2 and different values of gravity acceleration: 1 —g = 0 ; 2 —g = 
790 AV*/H = 3.232 x 10 11 cm/s 2 . Solid lines—Re 7 ; dashed lines— Im 7 . 


To illustrate this, in Figure 7.4, the dispersion curves are given in the cases of 
absence of gravity force (g = 0 ) and of large enough gravity force (g > g c ) when 
the aperiodic instability is suppressed. Within a wide range of evaporation 
front velocity V the value of stabilizing acceleration g c is nearly proportional 
to V 2 and, therefore, to 7 q, where Jo is the laser intensity. 

When k < mefi , the phase shift between the thermal and hydrodynamic 
perturbations (which manifests itself in the presence of the imaginary part of 
7 near the instability threshold) leads to the buildup of surface waves. At 
small evaporation rate, V mc/ix, this buildup may become stronger than 
the convective stabilization and evoke the instability. A characteristic form of 
the dispersion relation in the case of the coexistence of this and the aperiodic 
instability considered above is shown in Figure 7.5. Hydrodynamic oscillatory 
instability becomes important at relatively low radiation intensities when the 
short wavelength aperiodic instability cannot develop. At rather high gravity 
acceleration, the phase relation between thermal and hydrodynamic distur¬ 
bances is violated and instability is absent. 

Quantitative estimates of the growth rate and stabilizing gravity acceler¬ 
ation for hydrodynamic instability depend on the model used. In particular, 
a noticeable stabilization may be related to the melt viscosity that was dis¬ 
regarded in the previous analysis. This instability does not arise at a small 
film thickness either. It should be noted that a film thickness at which hy¬ 
drodynamic perturbations are stabilized is much larger than the radiation 
penetration depth or the thickness at which the melt flow becomes negligible, 
and the dispersion equation reduces to (7.21). In the present consideration, 
the description of the melt flow is simplified. We do not consider the melt 
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Figure 7.5: Dispersion curves 7 (k) at Jo = 2 x 10 7 W/cm 2 and different values 
of gravity acceleration: 1 —g = 0; 2 —g = 352 .6V*/H = 1.187 x 11 11 cm/s 2 ; 
3 —g = 1229 V^/H = 4.139 x ll 11 cm/s 2 . Solid lines—Re7; dashed 
lines— Im7. 


flow parallel to the melting front due to the recoil pressure gradient, and 
we consider the liquid layer thickness, H , as an independent parameter. In 
this situation, the description of hydrodynamic mode stabilization at a small 
melt thickness and the transition 7(if) —* 70 when H —> 0 has a qualitative 
character. 

Note that a finite thickness in the liquid film results in the reduction of 
the considered hydrodynamic instability and in the appearance of the new 
instability. The latter is connected to the nonpotential character of melt flow. 
The growth rate of this instability reaches its maximum at kH « 1. The dis¬ 
persion equation for this mode is similar to (7.22); however, the interaction 
of thermal and hydrodynamic perturbations leads to the growth of perturba¬ 
tion. The growth rate of this instability is usually small, as compared with the 
growth rates of the previous modes of instability. The acceleration required 
to stabilize this instability is also relatively small. 

As shown by calculations, a plane stationary evaporation front may be 
stabilized by a sufficiently strong gravity acceleration at arbitrary laser in¬ 
tensity. The dependence of the acceleration necessary to stabilize all of the 
instability modes considered above on the laser intensity is shown in Figure 
7.6. It should be noted that the relative position of the curves corresponding 
to different instability modes depends on the parameters of material, such as 
the surface tension, <t, and the absorption coefficient, /i. 

In typical cases of the coexistence of the short wavelength aperiodic and 
oscillatory hydrodynamic instabilities, the latter has a much smaller growth 
rate; but, to stabilize it, one requires a much higher gravity acceleration. A 
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/o. W/cm 2 

Figure 7.6: The stabilizing acceleration as a function of laser intensity. 

1— stabilization of short wavelength aperiodic (sublimation-type) instability; 

2— stabilization of oscillating hydrodynamic instability; 3—stabilization of 
long wavelength nonpotential hydrodynamic modes. 


stable evaporation regime can be reached in this case using a short interaction 
time (not sufficient for the development of relatively slow oscillatory instabil¬ 
ity) and the artificial gravity acceleration (to stabilize the aperiodic short 
wavelength instability). Note, however, that the case where the instability 
develops with short-pulse evaporation requires special study. 
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Chapter 8 

INSTABILITIES 
CONNECTED WITH 
LASER-INDUCED 
CHEMICAL REACTIONS 


In the previous chapters we considered several examples of instabilities appear¬ 
ing in macroscopic laser-matter interactions. In typical cases the mechanism 
of instability was connected with a positive feedback between the laser-induced 
temperature rise and the laser energy release. In its simplest form, this feed¬ 
back can be represented by an algebraic relation of the type of Equation (2.14). 
According to this equation, the absorption coefficient increases with the tem¬ 
perature increasing. This gives rise to the instability similar to the thermal 
explosion. 1 As we have seen while discussing the polymer ablation a feedback 
may have the form of differential equation (see Equation (2.20)). Further ex¬ 
amples of the feedback of this type can be found in laser thermochemistry. 2 In 
this chapter we shall consider two instability problems related to the surface 
oxidation of metals: (1) the instability of oxide film growth on the surface of 
uniformly irradiated target, and (2) the laser initiation of surface combustion 
waves. 


8.1 Instability of laser-stimulated 
surface oxidation of metals. 

Spatially uniform temperature fields 

In this section we shall consider the heterogeneous oxidation reaction of metals 
in air under the action of C02-laser radiation. We shall examine, first, the 
initial stage of heating of a metallic target, when the oxide layer is thin and one 
can neglect the interference phenomena leading to oscillations in the surface 
reflectivity. We shall assume for simplicity that the target is thermally thin, 
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i.e., that the characteristic time of temperature changes is much longer than 
H 2 /x, where H is the target thickness and x is the thermal diffusivity. Under 
this constraint, we may substitute the heat conduction equation with the 
equation of energy balance in the form: 

mcdT/dt = PA(h) - Pi{T) Pi = (3s(T - T 0 ) (8.1) 

where h is the thickness of the oxide layer, c is the specific heat of the oxide 
(per unit mass), P is the laser power, Pi is the power of convective energy 
losses, A(h) is the absorptivity of the layered system “metal -f oxide film”, m 
is the target mass, and s is the target area. We assume that the oxide film 
growth obeys the parabolic oxidation law: 3 

dh/dt = ( d/h ) exp (~T d /T) (8.2) 

where d and Td are constants. If h is small in comparison with the laser 
wavelength, the absorptivity A(h) can be written in the form (see Equations 
(8.8) and (8.9)): 

A{li) = Aq -f" bh 2 , (8*3) 

where Aq is the absorptivity of a metal. 4 The initial conditions for Equations 
(8.1) and (8.2) read 

T(0) = 0 and h( 0) = 0 (8.4) 

We introduce the dimensionless variables 

y = T/T d y 0 = T 0 /T d t = fist/me 
fi = 2bdPmc/T d (/3s) 2 y x = y 0 4- PA 0 //3sT d 

Under typical conditions of laser heating, the small parameters of the problem 
are y 0 1, yi <C 1, fi 1. Equations (8.1) and (8.2) can be reduced to a 
single differential equation 

d 2 y/dr 2 + dyjdr — /iexp(-l/y) = 0 (8.5) 

The initial conditions (8.4) are transformed into 

2 /( 0 ) = </o, dy/dt\ t =o = 2/i “ 2/o (8.6) 

The characteristic feature of the solutions of Equation (8.5) is an avalanche¬ 
like acceleration of the oxidation reaction after some induction period, 7*. 
This instability is similar to the inflammation process in combustion physics. 
The induction period usually is defined as the time at which the reaction rate 
reaches its maximum: d 2 h/dt 2 1 = 0. This time was calculated by Bunkin 

et al. 2 using the approximate solution of Equations (8.5) and (8.6) by the 
method of iteration. The result can be represented in the following form: 

n « (y x - y 0 ) -1 {1/ log[/x/(j/i - 2/o)] - 2/o} for g{n) » 1 
n « (j/i/m) exp(—1/2/1) forg(ri)<l 


(8.7) 
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A 0 ,% 


Figure 8 . 1 : Theoretical and experimental dependence of the induction time ti 
on metal absorptivity Aq for oxidation of copper in air by continuous 002 - 
laser. 1—the dependence calculated using Equation (8.7), g 1 (fast ac¬ 
tivation); 2 —the same for g 1 (quasi-stationary mode); the triangles— 
experimental results 5 . 


where g(r) = \/y(r) — \/y\. Thus, there are two modes of the oxidation 
reaction: fast activation, (g 1 ), and quasi-stationary activation (g 1 ). 
Figure 8.1 shows the comparison of experimental and theoretical dependences 
of the induction period on the initial absorptivity of the copper target. Ex¬ 
perimental data were taken from Arzuov et al ., 5 and the theoretical curves 
were calculated using Equation (8.7). 

A more general kinetic equation for the surface oxidation ( dh/dt oc h~ a ) 
was considered by Libenson . 6 The induction period there was calculated for 
the case of quasi-stationary activation. 

If the thickness of oxide film is comparable with the laser wavelength, 
interference variations in the absorptivity of the metal-oxide system can be 
observed. This is evident from an analysis of the equations for the absorptivity 
of a layered system 4 


A(h) = 1 — |r | 2 

1 — y/c 

ri2 = TT7^ 

Ao = 1 — |ri3 | 2 


r 12 exp (-2 iip) + r 23 
exp(-2iV>) + r 12 r 23 

1-v^ 

ri3 = -T, ~ 7- 

y/e — n — in 


</> = 2'Ky/e/X ) 


n 2 - n 3 

^23 = - 7 

n2H3 -1 

= n 0 - in , 



where r \2 and r 13 are the amplitude coefficients of reflection from the oxide 
and metal, respectively; e and €q are the dielectric permittivities of the oxide 
and metal, respectively; Ao is the absorptivity of the metal (Aq 1); h is 
the thickness of the oxide layer; and A is the laser wavelength. For thin oxide 
layers, A(h ) can be expressed as a series in powers of (h/ A): 
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A(h) = A 0 + a(h/X) + b(h/X) 2 + ... (8.9) 


where the coefficients a, 6 , etc. are determined by Equations ( 8 . 8 ). In the far 
infrared region we obtain 

a « 7Ti4o(n 2 — 1 ) and b « A'kA^tx 2 — 1 ). ( 8 . 10 ) 


Since Aq <C 1, it follows from (8.10) that a -C b. This means that the change 
in the absorptivity is determined by the quadratic term in the expansion of 
as assumed in (8.3). 

For weakly absorbing oxides (/iA < 1 , /i = 4ttk/X) a simple approximate 
expression can be derived 5 from ( 8 . 8 ) for A(h) 


A(h) 


n 2 A 0 -f 2n[ph - sin(p/i)] 
n 2 -f (1 — n 2 ) sin 2 (p/i/ 2 ) 


p = 47 rn/A. 


For the total number of interference oscillations, which can be observed in laser 
oxidation experiments, one can obtain a simple estimate (from the condition 
ph « 1 ): N ~ nj2nn. For example, for CU 2 O at the wavelength A = 10 . 6 /im, 
n = 2.45, k = 0.027, we have N « 14. A more complicated situation was 
analyzed by Bunkin et al ., 7 where the two-layer system (CU 2 O -f CuO) was 
considered. 


8.2 Instability of surface oxidation 

of metals: Two-dimensional effects 

In this section we shall show that the above considered one-dimensional growth 
of the oxide layer is unstable with respect to two-dimensional perturbations. 
We shall consider a thermally thin metallic plate and describe the temperature 
field by the heat conduction equation: 

dT/dt = xAT 4- A(h)I/Hcp (8.11) 

where x ls th e thermal diffusivity of metal, cp is the heat capacity (per unit 
volume), A is the Laplace operator in two dimensions (x,y), and A(h) is the 
absorptivity of the layered system (metal -f oxide). For simplicity, we will 
disregard the heat losses and assume the simplest linear oxidation law: 

dh/dt = dexp(—Td/T) (8.12) 

where d and Td are constants. If the laser intensity, /, is constant, Equations 
(8.11) and (8.12) have a particular solution T = T s (t), h = h s (t). We will 
study the stability of this solution with respect to small perturbations 

8T(y,t) = T\ exp 

8h(y,t) = h\exp 

Here, 7 (t) is a slow-varying function of time. A similar form of perturbations 
was used in Section 5.2, where the corrugation instability of nonstationary 


iky 4 / 7(r)dr 


0 


iky -f / 7 (r)dr 
0 


(8.13) 
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metal vaporization was studied. Calculations similar to those performed in 
Section 5.2 result in the following dispersion equation (see also references 8 
and 9). 

r . i 1/2 

7 (k) = -Xk 2 / 2 + [(xk 2 /2f + ( T dX I/KH)h dA/dh\ (8.14) 

where h = dh/dt. Due to the interference effects discussed in previous sec¬ 
tion, the derivative dA/dh is an oscillating function of h. One can see from 
(8.14) that the instability growth rate 7 (k) is positive for a positive dA/dh 
term. The maximum of 7 (fc) is reached at k = 0, i.e., the thermal conduction 
produces a stabilizing effect. As the oxide film thickness, h , increases, the 
derivative dA/dh changes its sign. When dA/dh is negative, the real part of 
the instability growth rate, Re 7 (fc), also becomes negative. This means that 
the perturbations, which were growing during the stage dA/dh > 0, begin 
to decay exponentially. We see, thus, that the stages of stable and unstable 
perturbation growth should change periodically. The evolution of monochro¬ 
matic initial perturbation is determined (in adiabatic approximations) by the 
value of the integral: 

r 

S(k,r) = / Re 7 (k,t)dt, (8.15) 

0 

where r is the laser pulse duration. The irreversible nonlinear stage of per¬ 
turbation growth is reached during the laser pulse, if the integral (8.15) is of 
the order of unity, S(k, r) > 1 . Consider the dispersion equation (8.14) in the 
case of low laser intensity, when 

(: T d xI/KHT 2 )h dA/dh <(xfc 2 /2 ) 2 (8.16) 

Expanding 7 (fc) in series in terms of small I and confining ourselves to the 
linear term, we have: 

7 (fc) « (T d I/ kHT 2 k 2 )h dA/dh (8.17) 

Substituting (8.17) into (8.15), we obtain 

5(fc, r) « AT d I/KHT 2 k 2 < 1 (8.18) 

The latter inequality can be proved in the following way. It follows from (8.12) 
and ( 8 . 11 ) that 

h/h » h/h w T d T/T 2 « T^AI/kHT 2 . 

Using this equation, the left side of Equation (8.16) can be rewritten in the 
form: 

(T dX I/KHT 2 )hdA/dh » {2nh/\){A 2 ){T dX I/ kHT 2 ) 2 (8.19) 

where n is the real part of the refraction index of oxide film; A is the wavelength 
of laser radiation; and (A 2 ) is the average value of A 2 (over several interference 
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Figure 8.2: The evolution of the oxide film on a thermally thin copper target. 
The values of time are indicated on the right. Normal incidence of CC>2-laser 
radiation, a —supercritical mode, b —subcritical mode. 


oscillations). The inequality (8.18) follows immediately from (8.16) and (8.19). 
We see, then, that at low laser intensities and arbitrary laser pulse lengths, 
the growth of the oxide film is stable, despite the fact that there are the stages 
when the perturbations (8.13) increase with time. 

We now will consider the case of high laser intensities, when 


(: T dX I/^HT 2 )h dA/dh >(xk 2 / 2) : 


( 8 . 20 ) 


If dA/dh > 0, the growth rate 7 (k) is real and increases as J 1 / 2 with increasing 
I. If dA/dh < 0 and the inequality (8.20) is fulfilled, 7 (k) is complex, and its 
real part, Re'y(k) = —xfc 2 /2, does not depend on I. It is clear that if the laser 
intensity exceeds some critical value /*, the amplification of perturbations 
during the stages with dA/dh > 0 should be greater than the damping of 
perturbations during the stages with dA/dh < 0. This means that the growth 
of the oxide layer is unstable. A rough estimate for J* follows from (8.20) and 
(8.19) : 8 

I m « {X/2nh) 1/2 KHT 2 k 2 /2(A)T d (8.21) 

A more precise determination of I* can be made using numerical calculations. 
Such calculations were performed by Gol’berg. 9 The oxidation kinetics were 
described by the parabolic oxidation law (8.2). This results in the dispersion 
equation 


l(k) = ~ (xk 2 + h/h)/2 

r . . 11/2 

+ y(xk 2 4- h/h) 2 /4 + {T d I X /KHT 2 )h dA/dh - X k 2 h/h 

which differs only slightly from (8.14). All the previously discussed qualitative 
results obtained using the linear oxidation law (8.12) remain valid for the 
parabolic law (8.2). In Figure 8.2 the nonlinear evolution of the shape of 
oxide film is shown. The pictures were obtained as the result of numerical 
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solution of Equations ( 8 . 11 ) and ( 8 . 2 ) for a thermally thin copper target. The 
initial conditions and parameters were chosen as: 

/i( 0 , x) = ho -f Sh ho = 1 x 10 “ 5 cm 6h = 0.25/io sin kx 
T( 0 , x) = 300 K I/H = 8.2 kW/cm 3 

case (a) k = 1.3 cm -1 and case (b) k = 5.3 cm -1 

Since the critical value of laser intensity depends on the perturbation wave¬ 
length (see ( 8 . 21 )), the same intensity I appears to be supercritical in case (a) 
and subcritical in case ( 6 ). Correspondingly, we see the rapid growth of the 
initial perturbation in the case (a) and the damping of the initial perturbation 
in the case (6). 

8.3 Laser-initiated surface combustion waves 

When laser radiation impinges on a target whose dimensions are greater than 
the focal spot size, heat transport effects play an important role. They can 
lead to a spatial instability in the form of a propagating reaction front. This 
situation is characteristic of heterogeneous exothermal reactions in which, un- 
ler certain conditions, a “hot spot” generated by laser irradiation begins to 
spontaneously expand. Let a laser beam with radius ro be incident on a 
metallic target occupying the half-space z > 0 . We assume that an exother¬ 
mal reaction with linear kinetics ( 8 . 12 ) occurs on the surface of the target. 
The temperature field in the target is described by the following boundary 
value problem : 2 

dT/dt = x [(1 /r)d/dr(rdT/dr) -f d 2 T/dz 2 ] 

—ndT/dz\ z=0 = AI(r) + pWdexp(-T d /T) 

T —> 0 for r, z —> oo 

Here W is the specific heat of the reaction, p is the density of the oxide, and 
d and T d are the constants in the linear oxidation law. We shall assume that 
the laser intensity distribution I(r) has the form 

I(r) = / TO exp(r 2 /rj)). (8.23) 

The analysis shows that the boundary value problem ( 8 . 22 ) and (8.23) has 
stationary solutions if the laser intensity I m does not exceed a certain thresh¬ 
old. The region in the plane (J m , ro) where the localized combustion regime is 
stable can be found from the stationary solution of Equations ( 8 . 22 ). Outside 
the region of stability, the problem ( 8 . 22 ) does not have stationary solutions. 

The stationary surface temperature near the axis of the beam can be writ¬ 
ten as: 

T(r,0) « T\(l — r 2 /R 2 ), r < R (8.24) 

The reaction energy release near the axis can be represented in a Gaussian 
form similar to (8.23): 

I c = pWdexp(—T d /T) « pWdexp(—T d /Ti) exp(—r 2 /b 2 ) 



(8.25) 
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Figure 8.3: The boundary of stability of localized oxidation in the plane (r, I). 
Above the boundary, a surface combustion wave appears. 


where b 2 = R 2 T\/Td . Since the heat conduction equation (8.22) is linear, its 
stationary solution can be written in the form: 

T(r,0) = To(r)+r c (r), 

where 

To (r) = (^/2)qT d exp(-r 2 /2rl)I 0 (r 2 /2rl) (8.26) 

is the temperature field arising from absorption of laser radiation (8.23), and 

T c (r) = (v / 7r/2)6/iTdexp(-T < i/Ti)exp(-r 2 /26 2 )/o(r 2 /26 2 ) (8.27) 

is the temperature field arising from the exothermic reaction (8.25). In (8.26) 
and (8.27), I${z) is the modified Bessel function of the zeroth order. We have 
introduced the parameters 

q = AI m ro/nT d and fi = pW dr 0 / kT<i 

Expanding (8.26) and (8.27) in powers of r near the beam axis and employing 
(8.24) and (8.25), we obtain a set of equations for T\ and R. It can be shown 
that this system of equations does not have real solutions for all values of q 
and pL. To find the boundary of the region where stationary solutions of (8.22) 
exist, we equate the determinant of the set of equations to zero. The result 
can be written 2 in parametric form as: 

q = 0.55(1 — a)(0.28 -f a 2 ) 

fi = 0.30(1 — a)(0.28 4- a 2 ) 3 (0.60 — a) exp 

where a is a parameter. Equations (8.27) determine the critical value of the 
dimensionless laser spot radius, f = ro(pWd/ nTd), as a function of the dimen¬ 
sionless laser intensity, I = Io(A/pWd). Figure 8.3 shows the boundary of 
the stability region. Below the stability boundary, the oxidation process is lo¬ 
calized spatially and the temperature field is stationary. Above the boundary, 
the high-temperature region expands with time. 


3.20(1 -fa) 

I (0.28 +a 2 )J 


(8.28) 
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Note that the appearance of a surface combustion wave has been observed 
experimentally by Arzuov et al., 10 where CC> 2 -laser stimulated combustion of 
tungsten in air was studied. At supercritical laser intensity (Jo ~ 50 kW/cm 2 ), 
a jump-like propagation of the combustion front was observed after the induc¬ 
tion period on the order of 1 s. 
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CONCLUSION 


The theoretical analysis summarized in this book shows that the processes of 
laser-matter interaction at moderate laser intensities can be unstable. The 
analysis of instability carried out previously was based on a simplified model 
description of laser absorption, energy transfer, and phase transitions in laser- 
heated material. It is important to understand to what extent the theoretical 
predictions agree with the vast amount of experimental material on laser- 
induced breakdown of dielectrics and laser ablation. 

The main focus in the previous analysis has been on the three different 
types of instabilities: 

(1) the hydrodynamic instability of vapor expansion into a vacuum or 
ambient gas; 

(2) the thermal instability (similar to the “thermal explosion”) arising due 
to the temperature dependence of the absorption coefficient; 

(3) the thermal instability of the phase transition front propagating in the 
direction of the temperature gradient. 

We now will discuss some of the experiments which can be interpreted as 
demonstrating the occurrence of these instabilities. 

Consider, first, the instabilities of vapor flow. When the vapor plume 
expands into an ambient gas, a shock wave is formed in front of the con¬ 
tact surface. In many cases of practical interest, the density of adiabatically 
expanding vapor is higher than that of the shock-compressed ambient gas. 
Under this condition, the contact boundary (whose velocity decreases with 
time) appears to be unstable. This instability, known as the Rayleigh-Taylor 
instability (RTI), leads to the distortion of the contact surface and turbulent 
mixing of the vapor and ambient gas. There are many experimental works in 
which this instability and mixing were observed. We would like to mention 
only the recent papers 1-4 devoted to excimer laser ablation of different materi¬ 
als (including high-T c superconductors). Spectroscopic and fast photographic 
methods were employed in these works. It was shown that at relatively high 
ambient gas pressure (po > 1 mbar) a strong mixing occurs at the contact 
surface. If the gaseous atmosphere is chemically active, the mixing stimulates 
the chemical reactions between the vapor and gas. In the case of high-T c 
superconductors the conditions are apparently well suited to the formation of 
diatomic metal oxides which are observed by spectroscopic methods. 3 Similar 
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effects appearing when detonation products expand into air were studied by 
Kull et al . 5 

Note that the Rayleigh-Taylor instability plays a very important role in 
processes of laser-driven inertial confinement fusion (ICF ). 7 Shell targets used 
in ICF experiments are susceptible to break up from RTI. The effects of RTI 
on the compression of ICF targets were studied extensively in the last two 
decades. However, the effects of RTI and turbulent mixing on the structure 
and optical properties of laser plumes were discussed mainly on qualitative 
level. More detailed experimental studies of these effects might be used to 
check semiempirical theories of turbulent mixing and to determine the values 
of empirical constants which enter into these theories. 

We turn now to the instability of thermal-explosion type. This instabil¬ 
ity arises if the volume energy release exceeds the energy loss due to heat 
conduction . 8 Maximum growth rate 7 (k) is reached at k = 0 when the heat 
conduction does not play role. There are two manifestations of this instabil¬ 
ity mechanism in laser experiments: ( 1 ) laser-induced thermal breakdown of 
transparent media, and ( 2 ) volume and surface thermochemical instability. 

A theoretical analysis of the breakdown of transparent media stimulated by 
small absorbing inclusions is described in Chapter 2 . We emphasize that the 
thermal-explosion behavior is connected with strong nonlinear temperature 
dependence of the absorption coefficient of dielectric. An inclusion plays the 
role of a seed, producing local heating of transparent material. Optical break¬ 
down of dielectrics containing small absorbing inclusions was studied in detail 
experimentally (references 9 and 10). Most of the experiments were conducted 
with optical glasses. The comparison of theoretical and experimental results 
conducted by Aleshin et al . 10 shows that calculated breakdown thresholds and 
induction periods are in reasonable agreement with experimental data. The 
model of inclusion-stimulated breakdown predicts the dependence of break¬ 
down threshold on the size of the focal volume. This dependence was observed 
experimentally both for gases and solid materials . 10-13 Since in many cases 
the characteristics of small absorbing inclusions are not well known, the in¬ 
verse problem arises: to determine the optical and statistical properties of 
the ensemble of inclusions by measuring the dependence of the breakdown 
threshold on the size of focal volume and on the laser pulse length. This 
problem was studied by Danilenko et al . 14 and Koldunov et al . 15 Note that, 
from a mathematical point of view, this problem is not well posed. Special 
mathematical methods have to be applied to solving problems of such type . 16 

Thermochemical instabilities are typical for C 02 -laser interaction with 
metals, semiconductors, and organic polymers. The laser-induced heating 
stimulates chemical reactions in the volume or on the surface of target. The 
rate of chemical reactions usually increases with temperature. If the products 
of a chemical reaction have a higher absorption coefficient than the initial 
substance, positive feedback arises, resulting in the instability. A well-known 
example of volume reaction leading to the thermochemical instability is the 
formation of soot as the result of thermal decomposition of organic poly¬ 
mers. Theoretical analysis of this instability mechanism was carried out by 
Tribel’skii and Liberman . 17 This analysis was developed further and experi- 
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mentally checked by Gol’berg et al. 18 It was shown that the model 17 provides 
a correct description of the instability development, and the calculated and 
measured values of the induction period are in good agreement. 

Laser-induced reactions of surface oxidation of metals were studied in great 
detail, both theoretically and experimentally. Interest in the study of these re¬ 
actions was stimulated by the technological applications of CC> 2 -lasers, whose 
radiation is highly absorbed by metal oxides, but is reflected from metal sur¬ 
faces. Various modes of oxidation and the influence of oxide layers on the 
dynamics of CC> 2 -laser heating of metals were examined (references 19-21). 
Variations in reflectivity during the laser-induced metal oxidation were used 
to study the kinetics of oxidation. 22 The instability of laser-induced oxidation 
of metals was first considered by Anisimov et al. 23 The results of numerical 
simulation of unstable oxide film growth taking into account the interference 
phenomena were described by Gol’berg and Tribel’skii. 24 The theoretical pre¬ 
diction of the thermo chemical instability was confirmed experimentally by 
Alimov et al. 25 The surface structures in the form of parallel bands whose 
period ranges from 13 to 200 /im, as well as localized disturbances of film 
thickness having a diameter about 100 /im were observed. A more detailed 
discussion of the problem of instability and structures formation, resulting 
from laser-stimulated heterogeneous chemical reactions, may be found in the 
literature. 26-28 

Now we will discuss the corrugation instability of a plane phase-transition 
front. This instability appears when the front of an energy-absorbing phase 
transition moves in the direction of temperature gradient. This instability is 
of thermal origin; it is expected to occur in various experiments connected 
with fast heating of condensed material (e.g., in shock wave compression of 
solids, fast Joule heating of metals, etc.). Similar instability is well known 
in the theory of crystal growth. 29 It was shown by Anisimov et al. 30 that 
plane vaporization front is unstable if the laser intensity exceeds some critical 
value. When the supercriticality is small, the nonlinear development of this 
instability leads to the formation of a stationary vaporization wave, moving 
at constant velocity. At higher intensities, the formation of nonstationary 
evaporation waves may be expected. Fluctuations of the recoil pressure as¬ 
sociated with this vaporization mode lead to the dispersal of target material. 
As a result, some part of the material is ablated in a condensed state. Many 
experimental facts are in qualitative agreement with this picture. The exper¬ 
imentally measured specific energy of laser ablation is several times less than 
the specific heat of vaporization. 31-33 An analysis of the phase composition of 
the ablation products 32 shows that a portion of the target material is removed 
in the liquid state. The measurements of target reflectivity (references 34-37) 
show that the total reflectivity, and especially its specular component, 37 is 
strongly reduced during the laser pulse. These facts can be explained by the 
instability of a plane vaporization wave. 

The above experiments may be considered as an indirect indication that 
the corrugation instability plays a role in the laser-matter interaction pro¬ 
cesses. As shown in Chapter 5, the spatial scale of perturbations having max¬ 
imum growth rate is of the order of light penetration depth. For metals, it 
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equals about 10“ 5 cm. Direct experimental study of this small-scale instabil¬ 
ity presents considerable difficulties. Anisimov et al. 38 studied the instability 
of thin film evaporation both theoretically and experimentally. In this case 
the spatial scale of growing perturbations is of the order of film thickness. 
Aluminum foils 10 to 50 /im thick were irradiated by CO 2 pulsed laser. Pe¬ 
riodic structures in the form of concentric rings were observed on both the 
front (irradiated) and rear sides of the target. The structures appear to be in¬ 
dependent of a change in laser polarization. Similar structures were observed 
by Zuev et al. 37 Complicated evaporation modes which may be interpreted as 
oscillatory and “spin” regimes were reported by Kuznetsov et al. 40 and Duley 
and Young. 41 

We have not considered structures of another type that were also detected 
on the surface of metals subjected to laser radiation. These structures re¬ 
semble a system of parallel bands whose period is on the order of the laser 
wavelength. The bands were oriented perpendicular to the direction of the 
electric vector of incident light wave (references 42 and 43). The most plausi¬ 
ble explanation of the origin of these structures is the interference between the 
incident electromagnetic wave and the surface electromagnetic wave generated 
at the metal-vacuum boundary. The theory of this phenomenon is reported 
by Bonch-Bruevich et al. 44 
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